VARTOOLS: A Program for Analyzing Astronomical Time-Series Data

This paper describes the VARTOOLS program, which is an open-source command-line utility, written in C, for analyzing astronomical time-series data, especially light curves. The program provides a general-purpose set of tools for processing light curves including signal identification, filtering, light curve manipulation, time conversions, and modeling and simulating light curves. Some of the routines implemented include the Generalized Lomb-Scargle periodogram, the Box-Least Squares transit search routine, the Analysis of Variance periodogram, the Discrete Fourier Transform including the CLEAN algorithm, the Weighted Wavelet Z-Transform, light curve arithmetic, linear and non-linear optimization of analytic functions including support for Markov Chain Monte Carlo analyses with non-trivial covariances, characterizing and/or simulating time-correlated noise, and the TFA and SYSREM filtering algorithms, among others. A mechanism is also provided for incorporating a user's own compiled processing routines into the program. VARTOOLS is designed especially for batch processing of light curves, including built-in support for parallel processing, making it useful for large time-domain surveys such as searches for transiting planets. Several examples are provided to illustrate the use of the program.


Introduction
Measuring the variation over time in the apparent brightness of a distant object is one of the primary techniques that an astronomer may use to study the universe. A few classic examples include using supernovae or pulsating variable stars as standard candles in setting the cosmic distance ladder (e.g. Hubble, 1929;Riess et al., 1998;Perlmutter et al., 1999), using eclipsing binaries to measure the fundamental physical properties of stars (e.g. Kopal, 1959;Torres et al., 2010), observing microlensing events to constrain the population of MAssive Compact Halo Objects (MACHOs) in the Galaxy or to identify extrasolar planets (e.g. Paczynski, 1986;Bond et al., 2004), and finding extrasolar planets by searching for stars that show periodic transit events (e.g. Henry et al., 1999;Charbonneau et al., 2000;Konacki et al., 2003).
Developments in technology, together with scientific interest in finding rare types of variability, have led to numerous massive time domain surveys, especially in the optical and infrared bandpasses (see for example Bakos et al., 2012, and references therein). Adding the time dimension to a survey leads to a massive flow of data, and a need for specialized software to analyze it. This need is expected to grow even larger as new time-domain surveys (e.g., LSST, TESS, PLATO, HATPI, etc.) continue to come online. This paper describes Vartools, a program which provides a suite of tools for processing light curves, and which is designed in particular for handling large data sets. An early version of Vartools, which included only a handful of period-finding tools, was briefly described in Hartman et al. (2008) who used it in analyzing the data from a variability survey of the open cluster M37. Since that time the program has been significantly altered, motivating the more detailed and up-to-date description provided here.
In addition to the numerous standalone software packages available which perform particular analysis tasks on light curves (e.g., transiting planet and eclipsing binary identification and modeling software; Kovács et al., 2002;Devor, 2005), there are several other software suites available for general astronomical light curve analysis. Notable examples include PyKE (Still and Barclay, 2012) 1 , the Peranso period analysis program (Paunzen and Vanmunster, 2016), the Period04 program (Lenz and Breger, 2005), Wqed which was developed for the Whole Earth Telescope project (Thompson and Mullally, 2009), and the open-source VStar program provided by the American Association of Variable Star Observers 2 .
The Vartools programs differs from these in a few key ways. First, it provides more general purpose analysis routines than the other programs. Second, unlike the other programs, it is designed specifically for processing large numbers of light curves. This includes the ability to seamlessly pass light curves between commands without requiring intermediate temporary files, not requiring any human interaction to process each light curve, and native support for parallel processing of light curves.
General signal processing packages exist for data analysis platforms such as matlab, R, idl, or python, but these are typically not tailored for astronomical purposes, and in general are collections of functions (like PyKE), requiring the user to solve a variety of bookkeeping problems, and in-so-doing to write a fair amount of code, to apply them to large astronomical time domain datasets. All of these platforms also allow executing system commands, so Vartools may be incorporated into pipelines built for these platforms as well.
The following section provides an overview of the program, including a general description of its operation, and discussions of the processing commands and control options that are provided; in section 3 several examples are provided; performance tests are presented in section 4; possible areas for future development are discussed in section 5. The appendices include comprehensive listings of the input and output syntax and data formats.

The Vartools Program
The basic operation of Vartools is to read-in a light curve, apply one or more "commands" (i.e. processing routines) to the light curve, and output resulting statistics to an ASCII table for later analysis. The program is written in ansi c, with regular updates provided on the web 3 . This article describes version 1.33 of the program, a static copy of which is preserved on the website. To allow maximum portability, it can be compiled with basic functionality without any non-standard external libraries, however a number of features have library dependencies (including the cfitsio library, Pence, 1999 4 ; the GNU scientific library, or GSL, Galassi et al., 2009 5 ; and the JPL NAIF CSPICE library, Acton, 1996 6 ). Compilation and installation are carried out using the GNU Build System 7 and has been successfully carried out on a number of system architectures (Linux, Mac and Windows). While Vartools may be used to analyze an individual light curve of interest, it has been developed primarily to be included in pipelines which conduct batch processing of large numbers of light curves.
As a simple example, listing 1 shows the use of Vartools to calculate the average, standard deviation, and expected standard deviation of the magnitudes in the light curve stored in the file "EXAMPLES/1" (this ASCII text file is included in the distribution of the program).
Listing 1: A simple example of running Vartools. In this example, and others shown in this paper, the text following "prompt>" is typed by the user into the shell. In all examples in this paper we assume a BASH environment. Note that the "\" character at the end of the line continues the command over multiple lines; the "oneline" argument, which causes one output statistic to be printed Alternatively, we can apply this to a list of light curves stored in "EXAMPLES/lc list" (each line in this list contains the name of a light curve file to process), and run this using 32 parallel processes (listing 2). A number of more complicated examples are provided in Section 3. The next subsection describes the general operation of the program, paying special attention to those features relevant for batch processing. The processing algorithms which are included as "commands" in Vartools are discussed in Section 2.2, built-in methods for extending Vartools with a user's own processing routines are discussed in Section 2.3, and options which may be used to control the operation of the program are discussed in Section 2.4. Performance benchmark tests are provided in Section 4.

Input
A user may provide either a single light curve or an ascii list of light curve files as input to the program. Each light curve is stored in a separate file, either in ascii text format (one measurement per row, with separate whitespace delimited columns for the time, magnitude or flux, uncertainties, and other optional information), or as a binary FITS table. The input list file may also be used to provide light-curve dependent data (e.g., this could be the right ascension and declination coordinates of each star) as needed by the various processing commands invoked by the user. There is flexibility for the user to specify the format of the input files. The user may also provide a string giving a set of commands that are passed to the shell and applied to the input light curve files before being read-in by the program (e.g., this may be used to process compressed light curve files without having to store the uncompressed files on the disk).

Data Flow
The input light curves are then sent to the commands in the order that they are invoked on the command-line. The user may choose from 51 built-in commands, as well as a potentially unlimited number of user-developed commands. In principle an arbitrary number of commands can be called in a single invocation of Vartools; though in practice a machine-dependent limitation will be set by the amount of available memory. These commands generally have parameters which control their properties (e.g., Listing 2: Processing a list of light curves in parallel. All quantities for a particular light curve are output on a single line, the "..." indicates that additional output from the program is not being shown here, and the parallel processing causes the output light curves to be in a different order from the input list. the range of frequencies to search in a periodogram, or the period, epoch, planet size, etc., for injecting a transit signal into a light curve). Most commands allow a number of options for specifying how the parameters are to be determined. Typically a parameter can be specified on the command-line for all light curves, it can be read-in from a column in the input list (so that a different value can be used for each light curve), it can be set using the results from a previously executed command (e.g., the period found by one command can be used in subtracting a harmonic signal with another command), it can be set to the value of a variable (see Section 2.1.4), or equal to the result from evaluating an analytic expression.
When batch processing hundreds of thousands (or more) of individual light curve files, often the speed of the process is limited by attempts to read-from or write-to the hard disk, rather than by computation speed. In such cases there is a significant advantage to eliminating redundant attempts to access the disk. This is one advantage of writing a pipeline in Vartools, where a light curve is read into memory once and output to the disk only when explicitly requested (output is treated as a command, which may be invoked at any time in the pipeline). When a similar pipeline is written using a shell script to piece together many standalone processing programs, it is often the case that the programs are not capable of transmitting the light curves via memory, forcing redundant disk access 8 . Moreover, scripting a pipeline around Vartools is generally simpler than using a variety of different programs which often require different data formats, necessitating messy, and potentially bug-prone data manipulation steps within a script. While pipelines written in an interpretted language such as matlab, idl, or python may also avoid redundant disk access, they generally have a slower execution than a compiled program like Vartools. Moreover building a pipeline from scratch in an interpreted language may well require more development effort for the user than using Vartools where all of the necessary bookkeeping is handled internally by the program. 8 in some cases this problem can be overcome by writing temporary files to shared memory rather than a physical disk.

Parallelization
One other significant feature of Vartools is its builtin support for parallel-processing of light curves. In the present version this is limited to shared-memory parallel processing, which is appropriate for processing on a single multi-core machine. The parallelization is done using the standard POSIX threads library, included with most C compilers. Vartools creates a number of threads specified by the user and assigns to each thread a light curve to fully process. Once that light curve is processed, the parent thread is notified, and a new light curve is assigned to the child thread. This is repeated until all light curves have been processed. The parallelization of individual commands (e.g. period searches) has not yet been implemented, but is under development. For the common case where the number of light curves exceeds the number of processors, the current parallelization implementation should be adequate. Note that when run in parallel mode, computed statistics are printed in the order that light curves finish processing. We use locks to prevent simultaneously writing out data for multiple light curves, with buffering of the output to minimize occurrences of multiple threads waiting for the same lock. As a result the ordering of light curves in the output table may not be the same as in the input list, and in general will differ each time the program is run.

Variables
Vartools supports the use of variables on the commandline to refer to vector or scalar data read-in from the light curves or list files, or computed by commands. These variables can also be used in evaluating analytic expressions with the -expr command, in controlling the flow of the program with the -if, -elif, -else and -fi commands, for processing multiple data columns within a single light curve with the -changevariable command, or in fitting analytic functions to the light curves with the -linfit or -nonlinfit commands. The special variable names id, t, mag, and err are reserved for a string giving the name of the image that a measurement comes from, the time, magnitude (or flux) and magnitude (flux) uncertainty, respectively. Additionally the special variable N R refers to the integer record number of a point in a light curve (1 for the first observation, 2 for the second, etc.). Finally a list of special names are reserved for analytic functions, constants, and operators (these are listed in Appendix C).

Commands
Here the built-in processing commands supported by Vartools are described. A detailed description of the expected input syntax is provided in Appendix A, while a description of the quantities output by each command is provided in Appendix B. The commands can be separated into several categories: (1) finding periodic signals; (2) calculating light curve statistics; (3) fitting models; (4) filtering; (5) simulating light curves; (6) manipulating light curves; and (7) controlling the flow of data through Vartools. In addition to the built-in commands, Vartools provides support for the dynamic inclusion of commands and functions developed by the user; this feature is described at the end of this section. Note that the choice of commands built-in to the currrent version of Vartools has been determined primarily by the authors' own research. As such there is a bias towards specialized routines that are useful in dealing with data from transiting planet surveys, while some general routines useful for other types of time-series applications may not be present.

Period-Finding Algorithms
Several of the commands included in Vartools are used to identify periodic signals in light curves. Each of these routines generates a spectrum giving the significance of a periodic signal as a function of frequency or period. The user can optionally output the spectrum of each light curve to a separate file, Vartools will also identify peaks in the spectrum and report the period or frequency, a measure of the significance, and other relevant information to the output ascii table. Below is a synopsis of each of these commands: 2.2.1.1. -LS. Calculates the Generalized Lomb-Scargle (L-S) periodogram of a light curve (Zechmeister and Kürster, 2009;Lomb, 1976;Scargle, 1982) using the method of Press and Rybicki (1989) to speed up the calculation of various trigonometric sums (we use the implementation of this "extirpolation" algorithm due to Press et al., 1992). The primary motivation for using the L-S method as opposed to directly applying the Fast Fourier Transform (FFT) to the light curve is that L-S handles nonuniformily sampled data (in astronomy this is the rule rather than the exception) without resorting to methods such as zero-padding or interpolation to a uniform grid, which amounts to making assumptions about the light curve at non-observed time instances. For this reason, this tool is widely used for identifying periodic signals in astronomy. The generalized periodogram extends the traditional L-S periodogram by allowing for a floating mean and non-uniform uncertainties (the user may optionally calculate the traditional periodogram instead; this is calculated following Press et al., 1992). The value of the periodogram at frequency f (in cycles per day) is given by equations 4 and 5 in Zechmeister and Kürster (2009), i.e.: where χ 2 0 is the value of χ 2 using the weighted mean of the data as the model, and χ 2 (f ) is the value of χ 2 using the best-fit sinusoidal signal with frequency f as the model. This statistic varies between 0 (for no signal present at all) and 1 (for a perfectly sinusoidal signal).
In addition to the peak periodogram value, the Vartools implementation of L-S outputs two measures of significance for peaks identified in the periodogram. These are the logarithm of the formal false alarm probability, and the signal-to-noise ratio (S/N).
There are several ways in which the false alarm probability may be estimated. We follow the Cumming et al. (1999) method described in section 3 of Zechmeister and Kürster (2009). In this method the signal is measured relative to the scatter in the residuals from the best-fit model (as opposed to using the scatter in the original light curve, or using the input measurement uncertainties, which are two other popular methods). Namely we take where N is the number of observations in the light curve, and LS best is the highest value of LS found in the periodogram. The false alarm probability is then given by where the false alarm probability for an individual trial frequency is and M is the number of independent frequencies sampled (sometimes called the "bandwidth" penalty). Here we differ slightly from Zechmeister and Kürster (2009) and use the Horne and Baliunas (1986) estimate of where f max is the maximum frequency calculated in the periodogram, and T is the time spanned by the light curve. Note that the false alarm probability will change if you scan a larger or smaller frequency range-this is expected behavior. It is also important to note that this false alarm estimate assumes Gaussian white noise. For real data the noise is often correlated in time, or does not follow a Gaussian distribution. In such cases the true false alarm probability is often much higher than the calculated probability, and should be determined through Monte Carlo simulations. Also note that M is a simple estimate that may be off by as much as a few dex in cases where the sampling is highly non-uniform. For applications requiring a false alarm probabilty accurate to better than 1 dex, it is necessary to conduct Monte Carlo simulations applying L-S to light curves with simulated noise to calibrate the FAP. The S/N measure provided by Vartools is determined within the periodogram as where LS is the value of the periodogram for a given peak, LS is the mean across the periodogram, and σ LS is the standard deviation of the periodogram. The user may optionally specify a clipping factor used to remove outlier points (peaks) from the values used in determiningLS and σ LS . The S/N value provides an alternative method for estimating the significance of a detection, however it is important to keep in mind that the "noise" in the periodogram does not follow a normal distribution, so the S/N itself does not follow a normal distribution (i.e., a spectroscopic S/N of 10 does not correspond to a 10σ detection). The user specifies the number of peaks to identify in the periodogram, these will be sorted based on the peakheight. If searching for multiple peaks, the user has the option to iteratively "whiten" the light curve after finding a peak (i.e., subtract the best-fit periodic sinusoidal signal from the light curve with the period fixed to the value identified in the periodogram), recompute the periodogram, and search for the next highest peak.

-GetLSAmpThresh.
Though not a method for identifying a periodic signal per se, this command is closely related to -LS and so it is discussed in this section. For many applications it is important to know not only that a signal is present in a light curve, but also what the minimum amplitude the signal could have had and still have been detected within the noisy observed data. This information helps in characterizing the selection effects present in a survey, and is thus vital for determining the completeness of a survey. This command fits a multi-harmonic Fourier series to a light curve at the period found by -LS. The model is then subtracted from the light curve, its amplitude is scaled by a factor α, and then it is re-added to the residual. The Lomb-Scargle FAP is calculated at the fixed period for this new light curve. The value of α is varied until the resulting FAP is equal to a threshold value chosen by the user; the resulting α times the observed amplitude is reported to the user.

-dftclean.
Calculates the Discrete Fourier Transform (DFT) of a light curve, and optionally applies the CLEAN deconvolution algorithm (Roberts et al., 1987). We follow the convention that the complex DFT of the set of points (t j , m j ), evaluated at frequency f is given by: and use the method of Kurtz (1985) to speed up the calculation. The user may also optionally output the window function (i.e. the DFT of the set of points (t i , 1) where t i are the times of observation, and the DFT is applied without subtracting the average from the y values). The DFT procedure takes two parameters: nbeam which determines the frequency sampling of the spectrum via ∆f = 1/(T × nbeam) where T is the time baseline of the observations, and the maximum frequency to calculate the DFT up to (the default is the Nyquist frequency given by 1/(2∆t min ) where ∆t min is the minimum time separation between consecutive points in the light curve). The CLEAN routine is controlled by the gain and the S/N threshold for selecting peaks in the DFT. After the DFT is calculated the spectrum is searched for the most significant peak above the S/N threshold. This peak, times the gain, is added to the cleaned spectrum, while the window function, scaled appropriately and shifted to the peak frequency, is subtracted from the DFT. The procedure iterates until no peaks above the S/N threshold remain.

-aov.
Calculates the Analysis of Variance (AoV) periodogram of a light curve, as defined by Schwarzenberg-Czerny (1989). The periodogram statistic follows a Fisher distribution of the form F (x; r − 1, n − r) where x is the value of the AoV periodogram, there are n data points in the light curve, and r bins are used in the fit. As with all period searches, one must also take into account the number of independent periods that are tried. We use the Horne and Baliunas (1986) estimate for the band-width penalty. Work requiring a more exact determination of the false alarm probability should conduct Monte Carlo simulations (see for example Paltani, 2004;Schwarzenberg-Czerny, 2012). Using the Fisher distribution to interpret the AoV statistic also naturally punishes overfitting the data. While the previous two commands use a sinusoid for the model signal, in this case the model which is fitted to the light curve is a discrete set of step functions (i.e., it is the phase-binned light curve). This method is thus comparable to the classical Phase Dispersion Minimization technique (Stellingwerf, 1978), but with a superior statistic. In principle -aov is more sensitive than -LS or -dftclean to light curves with sharp features, such as eclipsing binaries, for which a sinusoid is a poor approximation. However, because the phase-binned model is quite general (a relatively large number of free parameters must be used to fit typical signals), for most classes of variables the significance of the detection is lower when using -aov than another method such as -aov harm or -LS (e.g. Kovacs, 1980).

2.2.1.5.
-aov harm. Similar to -aov, this command computes an AoV periodogram for a light curve, however it uses a multi-harmonic Fourier series model for the signal, rather than the phase-binning model used by -aov. The fit is done using a projection onto orthogonal complex polynomials as described by Schwarzenberg-Czerny (1996.
Here one may also use the AoV statistic to optimize the number of harmonics to include in the fit.
2.2.1.6. -BLS. Computes a Box-Least Squares (BLS; Kovács et al., 2002) spectrum for a light curve. Unlike the other period-finding commands, -BLS does not search for signals that are continuously variable over the full phase curve, but instead searches for periodic box-shaped dips in a light curve. This algorithm provides an efficient means of identifying detached eclipsing binary light curves and transiting planet signals. The algorithm operates by scanning through frequencies, phasing the light curve at each frequency, and fitting to the phased light curve a model of the form for φ 0 + q ≤ 1, or for φ 0 + q > 1. Here φ is the phase,m is the average magnitude of the light curve, and δ m , φ 0 and q are free parameters representing the transit depth, starting phase of the transit, and transit duration in phase, respectively. The fit for the free parameters is done by binning the phased light curve into N bin bins, trying each of these bins for φ 0 and trying discrete values for q of the form i/N bin where i is an integer that varies between i min = floor(q min N bin ) and i max = ceil(q max N bin ), and q min and q max are set by the user. The user may optionally fit a simple "trapezoid" transit to the light curve for each of the significant periods found in the BLS spectrum. In this case the fractional duration of ingress and egress (i.e., in units of the total transit duration) is an additional free parameter in the fit. The trapezoid model is fit using the Downhill Simplex method (DHSX; Nelder and Mead, 1965) with initial parameter values for the frequency, transit depth, out-oftransit magnitude, transit time and transit duration taken from the BLS peak. The basic statistic computed by the BLS algorithm is the Signal Residue (SR) as a function of trial transit frequency, defined by Kovács et al. (2002) as where is the weighted sum of phased, weighted-mean-subtracted magnitudesx i within a trial transit (i.e., points with phase φ for a given frequency f , such that φ 0 < φ < φ 0 + q if φ 0 + q ≤ 1, or φ < φ 0 + q − 1 or φ > φ 0 otherwise),w i is the weight for point i normalized such that the sum of all weights in the light curve is unity, and we have r(φ 0 , q, f ) = transitw i . The maximum in equation 10 is taken over all pairs of φ 0 and q considered for a given transit frequency f . The function SR(f ) is also called the BLS spectrum. Kovács et al. (2002) also define the Signal Detection Efficiency (SDE) for a given peak in the BLS spectrum, given by whereSR is the average value of the BLS spectrum, and σ SR is its standard deviation. In many cases time-correlated noise leads to a BLS spectrum which slowly rises with decreasing frequency. To account for this, by default the Vartools program ranks peaks by a modified version of the SDE, called the spectroscopic signal to noise ratio (S/N), and calculated as follows. We first definẽ where the mean is computed, for a given frequency, over all trial combinations of φ 0 and q, and the quantity being averaged is the same quantity that is maximized in computing SR. A 3σ clipping is first performed before taking the mean. One may interpretSR(f ) as the average Signal Residue due to noise at a given frequency. The S/N is then taken to be where in this caseSR(f) is the average value ofSR computed over a frequency range from f −100∆f to f +100∆f with ∆f being the frequency sampling of the spectrum, and σS R is the standard deviation ofSR taken over the full spectrum. We perform an interative 3σ clipping before computing the average as a function of frequency and the standard deviation. This determination of S/N is similar to the SDE, except that it is the frequency-localized signal enhancement, relative to the global noise in the spectrum. Note that because SR does not follow a Gaussian distribution, this definition of S/N should not be interpreted using such a distribution. The user may optionally usē SR and σ SR in equation 14 rather thanSR(f) and σS R , in which case the S/N is the same as the SDE, except for the clipping applied in computing the average and standard deviation.
Another statistic which we calculate for each transit, and which may be used to select transit candidates, is the signal-to-pink noise ratio (S/N pink ) defined, following Pont et al. (2006), as: where σ w is the "white noise" r.m.s., σ r is the "red noise" r.m.s., n t is the number of points in transit, and N t is the number of distinct transits sampled. We take σ w to be equal to the r.m.s. of the light curve after subtracting the transit model, while we estimate σ r using the expression where σ bin is the r.m.s. of the residual light curve after binning in time with a bin-size equal to the duration of a transit, and σ bin,thy is the expected r.m.s. of the binned light curve if the noise were uncorrelated in time. In the event that the estimate for σ 2 r is less than zero, we take σ r = 0.
2.2.1.7. -BLSFixDurTc. Runs BLS on a light curve fixing the duration of the transit and the reference time of transit. Optionally the depth and the duration of ingress (for a trapezoidal-shaped transit) may be fixed as well. This may be used, for example, in a case where a single transit is observed with high precision (e.g., from space) and one wishes to search lower-precision data from another instrument for additional transits.
2.2.1.8. -BLSFixPer. Runs BLS on a light curve fixing the period (i.e., it finds the most transit-like signal in the light curve at a specified period). This may be useful for carrying out transit injection/recovery simulations to determine the detection efficiency where one does not wish to perform a full BLS search on the simulated light curves.
2.2.1.9. -wwz. Computes the weighted wavelet transform as defined by Foster (1996). This tool may be used to characterize signals with quasi-periodic behavior that evolve in time (in amplitude and/or frequency). The weighted wavelet transform is a discrete approximation to the continuous wavelet transform for the function x(t) given by for frequency ω, time-shift τ and wavelet kernel f (z). As shown by Foster (1996) the weighted wavelet transform is superior to the well known discrete wavelet transform for non-uniformly sampled time-series. Foster (1996) adopts the abbreviated Morlet wavelet (Goupillaud et al., 1984) for the kernel given by and approximates the integral using a weighted projection of the data onto the sinusoidal basis functions (e iω(t−τ ) ) with weights e −cω 2 (t−τ ) 2 . That is, the data are modeled as where the basis functions are and the coefficients y a are given by where the weighted inner product of functions f and g is given by with weights w α = e −cω 2 (tα−τ ) 2 (26) and the sums are over the observed times. The matrix S ab is given by Foster (1996) defines the Weighted Wavelet Transform (WWT) as which is a χ 2 statistic with 2 degrees of freedom with expected value 1. Foster (1996), however, suggests using the Z statistic (or Weighted-Wavelet Z-Transform, or WWZ) given by to characterize the significance of a signal with frequency ω and time-shift τ in the observed data. In both of these relations is the effective number of data-points contributing to the signal, and and are the weighted variations of the data and model. The amplitude of the signal (dubbed the Weighted Wavelet Amplitude, or WWA) is given by Given a value of c (by default c = 1/8π 2 is assumed), the -wwz command scans through frequency and timeshift (in keeping with the rest of the program we express frequency in cycles per day, whereas the frequency ω appearing in the equations above is in radians per day) computing the y a coefficients, WWT, WWZ, WWA and N eff . The user may optionally output this full transform (i.e., each of these quantities for every frequency and time-shift combination that is tested) to a file (as an ascii table, or a multi-extension fits image) for each light curve analyzed. For each trial time-shift the frequency with the maximum value of W W Z is also identified (i.e., f max,τ , and corresponding quantities W W Z max,τ , W W A max,τ , N eff,max,τ , and y a,max,τ ), and may be output to a separate file (this gives effectively the maximum signal as a function of time in the light curve). The maximum of W W Z max,τ across all trial τ values is determined, and its associated frequency, amplitude, N eff and model coefficients are included in the Vartools output ascii table. We also compute the median values of the max quantities across τ and include these in the output ascii table as well.

Light Curve Statistics
There are eight built-in commands which calculate statistics useful for charaterizing the variability of light curves.
2.2.2.1. -alarm. Calculates the statistic proposed by Tamuz et al. (2006) to detect variability. Letm be the weighted average magnitude of the light curve defined bȳ where the sum is over all points in the light curve. Then the alarm a is given by where a k is given by the sum on j is over the kth continuous set of points with the same value of sign(m −m), the sum on k is over all such values of a k , and the sum on i is over all points in the light curve.
and then binning this in lag time τ , such that the DACF at lag τ k is given by where the sum is over all ij pairs with lag τ ij = |t i − t j | satisfying |τ ij − τ k | < ∆τ /2, for a fixed sampling ∆τ , and M k is the number of these pairs. Equation 37 differs from the corresponding equation in Edelson and Krolik (1988) in that we use the product of the measurement uncertainties σ i σ j in the denominator, rather than (σ 2 − ǫ 2 ) with σ being the standard deviation of the time series, and ǫ being a "measurement error" associated with the data-set, and included to preserve the normalization. This modification allows for heteroscedastic data, but as a consequence, the autocorrelation at time-lag τ = 0 in general will not equal 1. An uncertainty on this correlation is estimated by taking the standard error on the mean for the points in a bin: (39) The DACF will be written out to a separate file for each light curve analyzed. This algorithm is sometimes used to identify variability periods for non-stationary quasiperiodic signals by searching for peaks in the DACF, however it is most useful in characterizing the coherence time scale of features (either true astrophysical signals, or correlated noise) in the light curves. Currently it is left up to the user to separately analyze the output autocorrelation function.

-chi2bin.
Calculates the reduced χ 2 after applying a moving mean filter to each light curve. This filter operates by replacing each m i in the light curve by the average value of m within a bin of size ∆t. That is, m i is replaced with where the sum is over points j such that |t j − t i | < ∆t. The light curve uncertainties after applying the filter are reduced by the amount expected for white noise. That is where σ ′ i is the new value for the uncertainty, and the sum is over the same points as above. One result of this is that if the light curve has time correlated variations, the reduced χ 2 reported by -chi2bin will increase as ∆t is increased. For pure white noise, the reduced χ 2 will be constant (to within statistical uncertainties) as ∆t is increased.
2.2.2.5. -Jstet. Calculates the J light curve variability statistic suggested by Stetson (1996), together with the Kurtosis K and the L statistic also defined by Stetson (1996). The J statistic is similar to -alarm in that its value will be greater for a light curve that shows timecorrelated variability than for a light curve with the same scatter, but for which the variations are not correlated in time. The J statistic is given by the expression where the sum in the numerator and denominator is of n pairs of observations, is the normalized residual for observation i, there are N observations total, observations i(k) and j(k) form pair k, and w k is a weight to be assigned to a pair. The Kurtosis of the light curve magnitudes (i.e., the fourth moment of the distribution) is estimated via and has a value of 2/π ≈ 0.798 when N is large and the magnitude values are drawn from a Gaussian distribution.
The L statistic is then defined as where w all is the maximum value of w that any light curve from the survey could have (i.e., a light curve with measurements from all images in a survey). This statistic provides enhanced significance over the J statistic for light curves with non-Gaussian magnitudes, and reduces the significance for light curves with missing observations. 2.2.2.6. -rms. Calculates the standard root-mean-square (RMS) scatter of the light curve given by: Additionally this command will compute the expected RMS of the light curve given the input photometric uncertainties. This is given bȳ where there are N p fixed frequencies f i in the fit, each having N harm,i higher-order harmonics and N subharm,i subharmonics, and the parameters m 0 , a i,k , b i,k , c i,k and d i,k are fitted by the procedure. Note that this same model can also be fit through the more general-purpose -linfit command; the separate -Killharm command is provided for convenience. By default the best-fit values for the parameters m 0 , a i,k , b i,k , c i,k and d i,k are reported in the output table. The user may optionally output amplitudes, phases, or relative amplitudes and phases instead. Here the amplitude of harmonic k is given by while the amplitude for sub-harmonic k is given by The amplitudes may also be specified relative to the amplitude of the fundamental frequency via and The phases are calculated as and where atan2(y, x) is the standard 4-quadrant inverse tangent function with output in radians. The relative phases of the harmonics and sub-harmonics are given by and Using relative amplitudes and phases for the harmonics and sub-harmonics is a common way to characterize the shape of a signal independently of the overall phase and amplitude. Such parameters are useful, for example, in classifying variable star light curves based on their morphology.
The -Killharm command also provides a measurement of the peak-to-peak amplitude of each harmonic series that is fit to the data (a separate amplitude is supplied for each period). To determine this amplitude, the model harmonic series is evaluated at 100×(N harm,i +N subharm,i + 1) phase steps. The minimum and maximum values of the series at the evaluated points are determined. Refined measurements are then performed near these points using the dbrent algorithm (Press et al., 1992).

-linfit.
Fits an analytic model that is linear in its parameters to a light curve. The user provides on the command-line an analytic expression together with a list of the names of variables in that expression that are to be varied. The fit is done by singular value decomposition (SVD) as implemented in the svdfit function of Numerical Recipes (Press et al., 1992). This function includes a correction which suppresses poorly constrained parameter combinations, with the goal of providing greater numerical stability. Because this can lead to differences from other tools which solve the linear least squares problem by direct inversion of the so-called design matrix, we provide a summary of the svdfit method below.
Let the light curve magnitudes be represented as a column vector m of length N , and suppose these are to be fit with a model mod which can be expressed as i.e., the linear combination of N param column vectors x i . These vectors are supplied by the user, and may be the times of observations (perhaps raised to some power), other columns read-in from the light curve (e.g., the X and/or Y coordinates of the star on each image), or any analytic expression involving such terms. The a i parameters are the N param variables to be fit by this command by minimizing χ 2 , which is given by Here the N × N param design matrix A has components with x ji being the ith component of vector x j and σ i is the ith measurement uncertainty in the light curve, and the vector b has components Applying SVD (via the svdcmp function in Numerical Recipes) to A yields an N × N param column-orthogonal matrix U, a N param × N param orthogonal matrix V, and a N param × N param diagonal matrix W with positive or zero elements, which together satisfy and which can be used to determine the least-squares solution for a via where U i and V i are the ith columns of each matrix, and w i is the ith diagonal element of W. The evaluation of equation 63 can be done efficiently by back-substitution (the svbksb function in Numerical Recipes). In cases where degeneracies exist between the different x i vectors, very small values of w i may be found leading to numerical instabilities in evaluating equation 63. This problem is overcome by setting 1 wi = 0 in cases where w i < TOL×w max , for some tolerance factor TOL and with w max being the maximum element in W. We adopt TOL = 10 −9 in Vartools.
Formal uncertainties on the parameters are also determined based on the covariance matrix-these uncertainties are accurate only if the measurement uncertainties supplied in the light curve are correct, and if the noise is uncorrelated in time and drawn from a normal distribution.

-MandelAgolTransit
. Fits the analytic model for the transit of a dark spherical planet across a limbdarkened spherical star described by Mandel and Agol (2002). This model is parameterized by the orbital period, an initial epoch of transit center, the ratio of the planet radius to the stellar radius, the ratio of the semi-major axis to the stellar radius, the normalized impact parameter (the minimum projected distance between the stellar and planet centers divided by the sum of the planet and stellar radii), the orbital eccentricity, the argument of periastron, and either two or four limb darkening coefficients (two for a quadratic limb-darkening law, and four for a non-linear limb-darkening law). Note that this choice of the impact parameter differs from some implementations where the normalization is done only to the stellar radius. Normalizing to the sum of the radii has the benefit of limiting the range of possible values to between 0 and 1 to ensure transit. The user may also simultaneously fit a Keplerian orbit to a set of radial velocity data, input as a separate file. In this case the orbital semi-amplitude and center-of-mass velocities are additional parameters. Because the transit model is nonlinear in its parameters, it is necessary to choose initial values from which to begin the fit. The DHSX is then used to search for a local χ 2 minimum near the initial position. The user may either fix the initial values on the command line, or use the results from the -BLS command to automatically initialize the parameters. Note that this routine is not intended to produce publication-quality fits, as it provides no mechanism for error determination, and does not support simultaneous fits of multiple light curves. It is instead useful for rapidly processesing many candidate transiting planet systems to identify initial parameter estimates (as might be done in a transit search pipeline).

-microlens.
Fits a simple microlensing model due to Wozniak et al. (2001) to a light curve. The model is given by where and The free parameters in the model include: the zero-point flux f 0 ; the lensed flux f 1 ; the impact parameter u 0 in units of the Einstein radius; the Einstein ring crossing time t 0 ; and the time of maximum lensing t max . The model is fit using DHSX.
2.2.3.6. -nonlinfit. Fits an analytic model to a light curve (non-analytic functions may also be fit through the use of user-developed functions, see Section 2.3.2). The model may be non-linear in its parameters. The user provides on the command-line an analytic expression together with a list of the names of variables in that expression that are to be varied, initial estimates for their best-fit values, and step-sizes to use in initializing the search algorithms. The user may optionally include a separate list of parameters that enter linearly into the model, and which will be optimized using linear least squares at each step in the nonlinear model fitting (via eq. 63). In addition to specifying the model and the free parameters, the user may also provide prior probability estimates for the parameters and/or constraints on the parameters (a list of analytic inequalities that must be satisfied by the parameters). The routine will attempt to minimize the expression: where m and mod are column vectors containing the measured and model magnitudes, respectively, Σ is the covariance matrix for the observational data, and P j are the N P priors which are analytic functions of the free parameters. By default Σ is taken to be a diagonal matrix with components Σ ii = σ 2 i , where σ i is the magnitude uncertainty for observation i. In this case the term ln(det(Σ)) is ignored by the routine as it is independent of the model parameters. The user may optionally specify an analytic expression for the uncertainties (this expression may include parameters which are varied in the fit), in which case the term will be included. Options for specifying a non-diagonal covariance matrix are discussed below. The value of X 2 is minimized over the region satisfying the parameter constraints. The quantity X 2 is related to the likelihood L via X 2 = −2 ln L + C for some constant C and is equal to χ 2 when using a diagonal covariance matrix and no priors. Note that the P j functions specified by the user are actually equal to −2 times the natural logarithm of functions proportional to the prior probability densities, the P j functions are not the prior probability density functions themselves. For example, if the user wishes to specify a Gaussian prior of the form: for some parameter a, they would express it as (a−ā) 2 /(σ 2 a ). We allow three different forms for describing non-diagonal elements in the covariance matrix (i.e., allowing for timecorrelated noise). These are the squared exponential, exponential and Matérn noise models. In the squared exponential model the covariance between times i and j is given by where σ i is the formal uncertainty for measurement i (if the user specifies an analytic expression to use for the uncertainty, then σ i is the expression evaluated at time i), δ ij is the Kronecker delta function, and a > 0, and ρ > 0 are parameters of the noise model. The exponential model is similar, with The Matérn model (Handcock and Stein, 1993) is given by where Γ is the standard gamma function, K ν is the modified Bessel function of the second kind, and a > 0, ρ > 0 and ν > 0 are free parameters. This model allows for covariances that fall off more slowly with separation than the squared exponential model. For ν = 1/2 the Matérn model is proportional to exp(−|t i − t j |/ρ), while for ν → ∞ it converges to the squared exponential. When non-diagonal covariances are used it is necessary to compute the inverse of the covariance matrix, which as an O(N 3 ) operation becomes computationally expensive when analyzing time series with tens to hundreds of thousands of points. To speed up the calculation we make use of the fact that for the noise models assumed we have Σ ij > Σ ik for j ≥ i and k > j if the light curve is sorted in time. Moreover, because the elements of Σ approach zero very rapidly when moving away from the diagonal, the covariance matrix is typically a sparse, near-diagonal matrix. We therefore compute the expression Σ −1 (m − mod) using Cholesky decomposition and back-substitution, truncating the loops within the algorithm when C ij < ǫ for an assumed tiny value of ǫ. The logarithm of the determinant of Σ is computed from this same decomposition.
To optimize the non-linear parameters the user may choose between using the DHSX method as implemented in the amoeba routine by Press et al. (1992), or running a Differential Evolution Markov Chain Monte Carlo (DEM-CMC) procedure.
The amoeba routine performs a greedy search for the X 2 minimum, and will tend to find the local minimum nearest to the initial starting position. The search is carried out using a set of N par + 1 vectors of parameters (the so-called simplex, here N par is the number of free parameters). In the Vartools implementation the components of each vector are initially set to the initial parameter values specified by the user, except that for vector j the component j − 1 is set to the initial value plus the specified step-size (for the first vector no component is adjusted). The routine evaluates X 2 for each vector of parameters, and then repeatedly adjusts the components of each vector following a prescribed method in an attempt to find the minimum value of X 2 . The routine stops once it achieves: where X 2 max is the maximum value of X 2 amongst the vectors in the simplex, X 2 min is the minimum value, ǫ = 10 −10 is a small number used to prevent division by zero, and TOL is the convergence tolerance (by default equal to 10 −10 , but the user may specify a different value on the command line).
The second option is to run a DEMCMC procedure (ter Braak, 2006). This procedure runs 2N par simultaneous Markov chains (a minimum of three chains are used). By default Vartools will first run the amoeba downhill simplex to find a local minimum about which the DEM-CMC will be executed, though the user may skip this step. One of the chains is then initialized to the best-fit parameter values, for the other chains we add to each parameter a Gaussian random number multiplied by the step-size. The initial user-supplied step sizes are iteratively adjusted to determine the step-size such that ∆χ 2 = 1 for each adjusted parameter, holding all other parameters fixed. We cycle through each of the chains proposing a new link N in the chain j using where k and l are the indices of two randomly selected chains, γ = 1.7/ 2N par , ǫ is a small number (0.1 by default, but can be modified by the user), r is a vector of Gaussian random numbers and U is an upper triangular matrix determined from the Cholesky decomposition of the parameter covariance matrix Σ (e.g. Press et al., 1992). We estimate Σ from the existing DEMCMC chains, recomputing U every 10000 links. We calculate ∆X 2 , the difference in X 2 between the proposed link and the prior link in the chain, and then accept the proposal with probability exp(−∆X 2 /2) (if ∆X 2 < 0, the proposed link is always accepted). If the proposed link is rejected, the previous link is copied to link N . The process terminates once either the total number of links reaches a user specified limit, or the total number of accepted links reaches a limit.
The N c chains are stored in a single array such that element i in this array is associated to chain i chain = i modulo N c . The user has the option of writing this array out to a file. By default all links are written, but the user may write out only a fraction of the links using the "printevery" option. To allow the user to interrupt the process without losing all of the results the links are written to the file buffer immediately after each MCMC step. However, we do not explicitely flush the file buffer, so the data will only be written to disk once the buffer is filled (this may be every few hundred links, but depends on the number of free parameters, and on the specific computer system used).
The full chain is also stored in memory and is used to compute statistics included in the output ascii table. By default these are the median and standard deviation of each free parameter. These may be interpretted as the most likely values and uncertainties for each parameter. The values corresponding to the lowest X 2 model are also included in the output ascii table. The user may change the statistics computed from the chain and may also specify different combinations of the parameters on which to compute the statistics. By default the first 10% of the chain is excluded from calculating the output statistics, however the user may optionally specify a different burnin fraction.
A limit is also provided on the total amount of memory allowed to be allocated to a particular MCMC command. Once this limit is reached, each time a new link is added to the chain the earliest link in the chain is forgotten. The default memory limit is 4 GB. Note that if Vartools is run in parallel mode, then the limit for each thread is the total limit divided by the number of threads.
2.2.3.7. -SoftenedTransit. Fits an approximate transit model due to Protopapas et al. (2005) to a light curve. The model is given by: where M 0 is the out of transit magnitude, δ is the depth of the transit, c is a constant controlling the "sharpness" of the transit, P is the period, T 0 is the reference time of mid transit, and η is the transit duration.
2.2.3.8. -Starspot. Fits a single, circular, uniform temperature spot model given by Dorren (1987) to a light curve. This model is given by: where M 0 is the out of transit magnitude, A and B are functions given by Dorren (1987) that depend on α, the spot angular radius, and β, the angle between the line of sight and the surface normal at the spot center, and a and b are quantities given by with µ ⋆ and µ s being the linear limb darkening coefficients for the unspotted photosphere of the star and for the spot, respectively, and F s /F ⋆ being the spot to photosphere flux ratio. Dorren (1987) suggest values of a = 0.0298 and b = 0.08745 for the Sun at 5000Å. The angle β is itself a function of the time of the observation, the rotation period of the star, the inclination of the stellar rotation axis, the latitude of the spot, and the longitude of the spot center at a reference time. This procedure is deprecated in the current version of Vartools, users should instead use the -macula model which is distributed with Vartools as a user-developed command.

Filtering
The built-in commands listed below can be used to filter light curves. Additionally, filtering may be done through the procedures described in the previous subsection (e.g., -Killharm to filter harmonic series from the data, or -linfit to apply a linear filter).

-clip.
Removes points from the light curve. Each point is removed from all vectors associated with the light curve (e.g. the time, magnitude, magnitude uncertainty, and other data read-in from the light curve). The default method is to σ-clip outlier magnitude values. In other words, given an average magnitudem measured RMS magnitude σ m , and specified clipping factor α, any point i with magnitude m i satisfying |m i −m| > ασ m will be removed from the light curve. The user may optionally run this in iterative mode, in which casem and σ m are both recalculated after clipping, and the clipping is repeated with the new values; this can either be repeated for a specified number of iterations, or until no further points are removed. Finally the user may also simply remove points with formal errors ≤ 0 or undefined magnitude values. This later capability can be used in combination with an -expr command to perform a customized clipping procedure (e.g., listing 3).

-medianfilter.
For each point i, determines the median magnitude value of all points j that satisfy |t i − t j | < ∆t for fixed ∆t, and then either subtracts this value from the magnitude of point i (high-pass filtering the light curve), or replaces the magnitude of point i with this value (low-pass filtering the light curve). The user may optionally use the average, or the uncertainty-weighted average instead of the median.

-SYSREM.
Performs the trend-filtering procedure suggested by Tamuz et al. (2005), which effectively amounts to determining the principle components of the light curve dataset, and filtering these from the light curves. The filter operates by attempting to model the light curves indexed by i consisting of points indexed by j as: where ǫ i,j is a gaussian random number drawn from a normal distribution with meanm i and standard deviation σ i,j ; and there are N trend trend vectors a k,j , called "airmass-like" terms by Tamuz et al. (2005), to which each light curve has a linear response parameterized by c k,i , called "color-like" terms by Tamuz et al. (2005). Initial airmass and/or color-like vectors are assumed (these might be the airmasses of the images, or the x/y pixel positions of the stars, for example). For each assumed airmass-like (color-like) vector the associated color-like (airmass-like) vector is solved for by linear least squares minimization. The so-determined color-like (airmass-like) vector is then fixed, and the airmass-like (color-like) vector is then redetermined by linear least squares minimization. The procedure is repeated for the next vector until all N trend trends are fitted. To prevent outliers from dominating the determination of trends, the user specifies two σ-clipping parameters. The first for excluding points from individual light curves in determining the average magnitudes of the light curves, the second for excluding points from the fits for the airmass-like or color-like terms.
2.2.4.5. -TFA. Applies the Trend Filtering Algorithm suggested by Kovács et al. (2005). This is similar to SYSREM in that the light curves are modeled as a linear combination of trend vectors as described in equation 80, however in this case the a k,j "airmass-like" terms are taken to be N trend template light curves chosen randomly from among the full ensemble of observed light curves, and only the c k,i "color-like" terms are determined by fitting. One other practical difference between SYSREM and TFA is that the latter typically uses more trend vectors in modelling the light curves (e.g. for a field containing 10000 stars, each containing ∼ 10000 observations, it is not uncommon to use ∼ 1000 template light curves in TFA, whereas typically applications of SYSREM may only use a few dozen trend vectors). When the same set of templates is used to filter numerous light curves the filter may be applied efficiently by a combination of Singular Value Decomposition (or alternatively Lower-Upper Diagonal Decomposition) applied once to the matrix of templates, followed by the much faster Singular Value Back Substitution procedure applied to each light curve to be filtered.
2.2.4.6. -TFA SR. Applies the Trend Filtering Algorithm in "signal-reconstruction" mode. In this case the assumed light curve model is modified from that given in equation 80 to: where S i (t j ) is a light-curve-dependent signal. In Vartools the allowed signals include: (1) a periodic stepfunction of the form where Θ(X) is the Heavside-step function (= 0 for x < 0; = 1 for x > 0), P is the fixed period of the signal, there are N bin steps used, and the free parameters are the α u terms; (2) a fixed signal form evaluated at the times of observation (and read-in from a file); or (3) a harmonic series of the form given in equation 49. Additionally the user may include a number of light-curve-specific vectors among the series of trend vectors used in the fit. In practice this slows down the filter as it requires a re-inversion of the design matrix used in the linear fit.

Light Curve Simulation
The following built-in commands may be used in simulating light curves.
2.2.5.1. -addnoise. Adds noise with optional time-correlation to the light curve. Five different methods may be selected from. The first is simple uncorrelated Gaussian white noise, in which case the user specifies the standard deviation of the noise.
The second, third and fourth methods simulate a Gaussian process with a fixed covariance matrix Σ. These methods differ in the form of the covariance matrix. The second method assumes a squared exponential model where the elements of the covariance matrix are given by Eq. 69. The third method assumes an exponential covariance with elements given by Eq. 70. The fourth method uses the Matérn class of covariance matrices with elements given by Eq. 71. In all three cases the Cholesky decomposition of Σ is computed and then multiplied against a vector of random numbers drawn from a normal distribution to produce a light curve drawn from the specified covariance matrix.
In the fifth method the noise is modeled as where ǫ w is a white-noise component that is a random number drawn from a Gaussian distribution with zero mean and standard deviation σ w (specified by the user), while ǫ r is a red-noise component drawn from a distribution with standard deviation σ r and time-correlations such that the power-spectral density of the noise scales as 1/f γ (both γ, which must satisfy −1 < γ < 1, and σ r are specified by the user; taking γ = 0 yields uncorrelated noise). We use the method of McCoy and Walden (1996) to generate the red-noise. The McCoy and Walden (1996) procedure for simulating a series of N = 2 p points uniformly sampled in time, described in Section 4.2 of that paper, works as follows: (1) generate a set of points d m,k with m = 1, . . . , p and k = 1, . . . , 2 p−m each drawn from a Gaussian distribution with zero mean and variance S m , and generate an additional point x p , 1 with zero mean and variance S p+1 (McCoy and Walden, 1996 gives a prescription for determining the variances S m that depends on σ r and γ); (2) apply the inverse discrete wavelet transform (we use the implementation in the GNU scientific library) to the vector (x p,1 , d p,1 , d p−1,1 , d p−1,2 , . . . , d 1,1 , . . . , d 1,2 p−1 ). To generate non-uniformly sampled observations, we first simulate a series of N = 2 p uniformly sampled points as above, taking p = min(20, ceil(lg(T /δ t ))) where T is the time spanned by the non-uniformly sampled observations, and δ t is the minimum time between successive observations. We then linearly interpolate this sequence onto the non-uniformly sampled observation times.
2.2.5.2. -copylc. Makes copies of the light curve under analysis. Data read in from the input list file, as well as data computed by commands executed prior to the copy command are also replicated. Each copy is then independently processed through the subsequent commands issued to Vartools. This command is useful, for example, for carrying out simulations in which different signals are injected into the same light curve to estimate the recovery efficiency, or different pure-noise light curves with the same time sampling are simulated to determine the bandwidth correction for period finding algorithms such as -LS.

-Injectharm.
Adds to the light curve a harmonic series of the form where P is the period of the signal, A 1 , A k , and B k are the amplitudes of the fundamental, harmonic, and subharmonic terms, respectively, and φ 1 , φ k , and ψ k are the phases of the fundamental, harmonic, and subharmonic terms, respectively. The user indicates the number of desired harmonic and subharmonic terms (N harm and N subharm ) and how each of these are to be determined (their values may be fixed on the command-line, read-in from the input list of light curves, or randomly drawn from uniform or uniform-log distributions). The amplitudes and phases of the harmonic and subharmonic terms may also be specified relative to the amplitude and phase of the fundamental (i.e. the input is R i1 = A i /A 1 rather than A i , and/or φ k1 = φ k − kφ 1 rather than φ k ). This latter option is useful if one wants to inject into the light curve a signal with a fixed shape (e.g. a saw-tooth like pulsation), but random overall amplitude and phase, which someone might want to do to determine the completeness of a variability survey to a particular class of variable stars. In this case one would fix the relative phases and amplitudes of the harmonic and/or subharmonic terms, but allow the amplitude and phase of the fundamental to vary.  Mandel and Agol (2002) model for the transit of a nonluminous spherical object in front of a limb-darkened spherical star. The user may specify the model parameters, including: the orbital period, the radius of the planet, the mass of the planet, the phase of the transit at time T = 0, sin i (i is the orbital inclination), the eccentricity and argument of periastron or e cos ω and e sin ω, the mass of the star, the radius of the star, quadratic or non-linear limb darkening coefficients for the star, and an optional light dilution factor. Alternatively some of these parameters may be drawn from random distributions. For the period this can be a uniform distribution in P , log P , f = 1/P or log f . For the planet mass and/or radius this can be a uniform distribution or a log-uniform distribution. For sin i it is a uniform distribution of orbit orientations in space, subject to the constraint that there must be a transit, which corresponds to cos i being drawn from a uniform distribution between 0 and C with

Light Curve Manipulation
The following commands can be used to perform common light curve manipulation tasks: 2.2.6.1. -binlc. Bins the light curve in time. The user either specifies the number of bins to use, or the size of the bins in units of the time coordinate, and may also control the starting time of the first bin (by default it is the time of the first observation in the light curve). The user may either use the average, median or error-weighted average for the resulting binned value. All vectors associated with this light curve (i.e. columns read-in from the light curve file including the magnitudes), except for the time and uncertainty, vectors storing string or character data, or any vectors explicitly indicated by the user, are binned as follows: while the uncertainties are binned following: wherex i is the ith binned value of vector x, andσ i is the ith binned uncertainty, the sum on j in each case is over the unbinned points with times that fall within bin i, N i is the number of points that contribute to bin i, and the uncertainties are the standard errors on the binned quantities assuming Gaussian uncorrelated noise. Bins for which N i = 0 are excluded from the binned light curve. The time reported for each bin is either the time at the center of the bin, the average of the times that fall within the bin (i.e.t i = j t j /N i ), the median of the times that fall within the bin, or the user may choose not to shrink the size of the light curve, in which case all points within a bin are replaced by their binned value, but the light curve times are not changed. The user may also optionally provide a list of variables which will not be binned in the default manner. For each variable in the list the user also indicates which statistic to use in calculating the binned value. Options are the same as for the -stats command.
2.2.6.2. -changeerror. Sets the magnitude uncertainties σ j in a light curve equal to the measured r.m.s. scatter of the light curve magnitudes.
2.2.6.3. -converttime. Converts the time-system of the light curves between julian date (JD), modified julian date (MJD=JD−2400000.5), helio-centric julian date (HJD; i.e., JD approximately corrected to the center of the solar system assuming observations are made at the Earth-Moon barycenter, and that this position follows an elliptical orbit about the center of the solar system with linear perturbations to the orbital elements), and bary-centric julian date (BJD; i.e., time corrected to the center of the solar system using the NASA JPL ephemeris to determine the location of the observer relative to the barycenter at a given time); constants may be added to or subtracted from the input and/or output times. Time conversions between a UTC-based system (i.e., JD values having been calculated directly from UTC without accounting for leap-seconds; this is not formally correct, but a very common practice) and a TDB-based system (barycentric dynamical time, i.e. the UTC time is corrected to a leap-second-free time system before converting to JD) are also supported. In the year 2015 there is approximately a one minute difference between the two systems (see the discussion in Eastman et al., 2010).
Conversions to/from BJD and between UTC and TDB use the NASA JPL NAIF cspice library (Acton, 1996) following the expressions given by Eastman et al. (2010); we do not allow for Shapiro or Einstein time-delays in converting to BJD (i.e., gravitational redshifts are ignored).
Differences between BJD and HJD are due mostly to the orbits of Jupiter and Saturn (which are ignored for HJD), and can be as large as 4.2 seconds for observations made between 1980 and 2020.
The user may specify an observatory where the data were collected, or the latitude, longitude, and altitude. If this is not done then the center of the Earth is assumed, which can introduce up to a 21 millisecond error in the BJD correction. The procedure optionally accounts for proper motion and/or coordinate precession.
Based on comparing time conversions made with Vartools to those computed through the JPL HORIZONS interface 9 we conclude that this command has an internal precision of ∼ 0.1 milliseconds for conversions near J2000.0. At present this command does not support time conversions for observations made off of the surface of the Earth.
2.2.6.4. -difffluxtomag. Converts the light curves from differential flux units (which are output, for example, by the popular ISIS image subtraction procedure; Alard and Lupton, 1998) into magnitudes. The expression for converting from differential flux df to magnitude m is given by: where m 0 is the magnitude of a source yielding a flux of 1, ∆ m is an offset that the user may give to switch to a different magnitude system, and the reference flux f ref is determined from the input reference magnitude m ref via: 2.2.6.5. -ensemblerescalesig. Transforms the magnitude uncertainties of the light curves by the expression Where σ ′ i,j is the new uncertainty for light curve i, observation j, after the transformation. Here a and b are determined by fitting a linear relation of the form: where for light curve i, and for light curve i, andRM S is the expected RMS of the light curve (equation 48). The result of this transformation is that, for many datasets, the χ 2 i /dof values are distributed about unity. The justification for this is given in Appendix D. The parameter a is related to the ratio of the gain and the effective gain, while b is related to a constant magnitude error-term (e.g. Hartman et al., 2005). This routine requires all of the light curves to be read-in simultaneously, and thus is incompatible with the -parallel option in the current version of Vartools.
2.2.6.6. -expr. Sets a variable equal to an analytic expression. For example, one might give the command "expr mag=mag/2" which evaluates mag/2 at each point in the light curve and replaces the magnitude with the result. If the variable on the left-hand-side of the equality has not previously been defined, it will be created. Variables which appear on the right-hand-side can be the name of a vector which is read-in from the light curve, a scalar or vector created by another command (e.g. the fitting parameters, or the output model vector, created by the -linfit command), or any output parameter from a previously executed command (such as the light curve RMS computed with the -rms command). Vartools uses an implementation of the precedence climbing algorithm (Richards and Whitby-Strevens, 1979) for parsing and evaluating analytic expressions. The list of recognized functions is provided in Appendix C, while a method for users to define their own functions which may then be used in such expressions is discussed in Section 2.3.2.
2.2.6.7. -fluxtomag. Converts a light curve from fluxes into magnitudes. The uncertainties are also converted from flux uncertainties into magnitude uncertainties.
2.2.6.8. -Phase. Replaces the time coordinate of a light curve with its phase (from 0 to 1) and sorts the light curve by the phase. The user must specify the period to use for phasing the light curve, this may either be fixed on the command line for all light curves, read-in from the input light curve list, or taken from a previous command (e.g. the period found by -LS). The user may also optionally control the reference time for phase 0.
2.2.6.9. -resample. Resamples the light curve onto a new time base. Several methods are allowed for interpolating or extrapolating the data, including setting resampled points to the value of the observation that is closest in time, performing linear interpolation between points, cubic spline interpolation (e.g., Press et al., 1992), cubic spline interpolation with a constraint that the interpolating function be monotonic between points (Steffen, 1990), and Basis-spline interpolation (implemented using GSL routines). The user may optionally adopt different interpolation methods for resampled points that are close to observations and those that are far from observations. 2.2.6.10. -rescalesig. Rescales the magnitude uncertainties for a light curve (i.e. σ = ασ) such that χ 2 /dof = 1 for that light curve.

Control of Data Flow
The following commands control the flow of data through the pipeline: Changes the variable used for the time, magnitude, magnitude uncertainty, or image identifier in subsequent commands. For example, the command "-changevariable mag mag2" would cause subsequent commands to Vartools to operate on the data stored in the variable mag2 for the light curve magnitudes. The data stored in the original variable mag remains unchanged, and may be used in analytic expressions.

2.2.7.2.
-if, -elif, -else, -fi. Can be used to make execution of commands conditional upon the evaluation of expressions, following the typical logic of "if", "else if", "else" and "end if"-type conditional statements. The use of the terms -elif and -fi are as in the Bourne-Again SHell. Vartools supports nested conditionals, but all conditional constructs are ignored by commands which process all light curves simultaneously (e.g. -SYSREM) as well as by the -savelc and -restorelc commands.
Outputs the light curve in its present state to a file (or to stdout if specified by the user). The user can specify a rule for naming the output file, and can also control the format of the output data. By default the light curve will be output as an ascii text file, but may also be output as a binary FITS table.
2.2.7.4. -restorelc and -savelc. Are used to save a light curve (including all associated vectors) and later restore it to this previously saved state.

Miscellaneous Commands
The commands listed here do not fall under the other categories.
2.2.8.1. -findblends. This command may be used to determine if a signal detected in a given light curve is likely due to blending with a nearby variable source. The command follows the procedure described by Hartman et al. (2011). Each light curve in the input list is spatially matched to another list of light curves (by default the input light curve list is matched to itself, but the user may optionally provide a different list for the matching). A harmonic series with a specified period (either read-in from the input list file, fixed on the command line, or set equal to the output from a previously executed command) is then fit to both the input light curve, and to all matching light curves. The peak-to-peak amplitude of the fitted harmonic series is then determined for each source and converted to flux units (by default the routine assumes the input is in magnitudes). The light curve with the highest flux amplitude is reported. If this is not the same as the input light curve, then the neighboring object is the most likely source of the variability.

User-Developed Extensions
While Vartools includes a wide-range of built-in processing routines, it falls far short of encompassing the full range of routines that one might wish to apply to a light curve. It is therefore essential to enable a user to incorporate his/her own light curve processing algorithms within the program. To that end two mechanisms have been included: one for dynamically loading at run-time a compiled library defining a user-developed command; and another for dynamically loading a library defining a userdeveloped function which can be used in analytic expressions on the command-line. Each of these are described in turn.

User-Developed Commands
Each library implementing a user-developed command must contain five functions, with standardized names and calling syntaxes, which are used to interface the command with Vartools. These functions perform such tasks as initializing variables defining standard properties of the command, parsing the command-line when the user calls the command, displaying the syntax and help information for the command, and executing the command on a light curve. These required functions should be written in C, while the algorithm itself may be written in any language callable from C. Templates are provided of the functions required in each library, together with a Makefile for compiling and linking the code into a dynamically loadable shared object library. A template is also provided illustrating how to include a processing algorithm written in F ORT RAN in Vartools.
Vartools uses the GNU Libtool package 10 for compiling and dynamically loading libraries in a system-independent fashion. When Vartools encounters an unrecognized term on the command-line, it will search for an appropriate library in a data directory defined during the installation of Vartools, load it if available, and execute the command. So long as the libraries are installed in the correct directory, they may be called by the user in the same way that all other commands are executed. To avoid performing a directory search, the user may also use the -L option to explicitly load a library file.
A number of example libraries are distributed with Vartools which provide additional commands. These include: 2.3.1.1. -fastchi2. Calculate the Fast χ 2 periodogram following Palmer (2009). This routine scans through a range of frequencies fitting a harmonic series to the data at each trial frequency and reporting the minimum χ 2 . It uses the fast Fourier transform to do the fit quickly in O(N log N ) time. The user specifies the number of harmonics to include in the series, the minimum and maximum frequencies to search, the order and reference time of a polynomial in time which is removed from the light curve before performing the search, the time-span to use in computing the Nyquist frequency (if not calculating it automatically), an over-sampling factor for determining the resolution of the periodogram, and a tolerance parameter used in performing a higher resolution search around the peaks in the periodogram.

-jktebop.
Fit or inject a detached eclipsing binary light curve model [in]to a light curve. This command uses the code from the JKTEBOP program due to Southworth et al. (2004), converted into a callable library. The JKTEBOP code is in turn based on the Eclipsing Binary Orbit Program (Popper and Etzel, 1981;Etzel, 1981) and implements the model of Nelson and Davis (1972). The JKTEBOP-based library is written in F ORT RAN , so this command can be used as a template for incorporating a F ORT RAN algorithm into Vartools. The basic parameters for the model include the orbital period of the binary, a reference primary eclipse epoch, the normalized sum of the component radii ((R 1 + R 2 )/a), the ratio of the radii (R 2 /R 1 ), the mass ratio (M 2 /M 1 ), the surface brightness ratio (J 2 /J 1 ), the orbital inclination or normalized impact parameter, the eccentricity parameters e sin ω and e cos ω, limb darkening coeffections, gravity darkening coeffections, reflection effect coeffections, third light, and a tidal lead/lag angle. The fit is performed using DHSX.
2.3.1.3. -macula. Fit or inject a Macula spot model [in]to a light curve. The model and F ORT RAN code implementing it are both due to Kipping (2012). This model allows for an arbitrary number of starspots which are allowed to evolve over time with a linear growth/decay law. The parameters in the model include: the equatorial rotation period; the inclination of the star; quadratic and quartic differential rotation coefficients; four limb-darkening coefficients for the stellar photosphere and four for the spots; and a blend parameter. Additionally each spot is parameterized by: its longitude at the time of maximum size; its latitude at the time of maximum size; its maximum angular size; the spot-to-star flux ratio; the time of maximum size; the lifetime; the duration of ingress (time to appear); and the duration of egress (time to disappear). The model can be fit either with DHSX, or with the Levenberg-Marquardt algorithm (Press et al., 1992). The latter makes use of partial derivatives that are provided by Kipping (2012).

-magadd.
Adds a constant to a light curve. This is provided as a trivial example of a user-developed command, and may be used as a template.

User-Developed Functions
The mechanism for incorporating a user's own functions into Vartools is similar to the mechanism for incorporating user-developed commands. Again the user creates a dynamically loadable library written in C. The library must contain an initialization function with a standardized name which is called when the library is opened. This function is used to indicate the names of all new analytic functions supplied by the library, the number of double-precision arguments expected by each function (at present all functions input and output double-precision variables), and pointers to the C functions in the library (i.e., variables storing their location in memory) which execute the computations associated with each analytic function. These C functions are expected to take a vector of double-precision numbers, with each component in the vector representing a separate argument to the function, and to output a single double-precision quantity as the result. The functions may in turn be wrappers to routines written in other languages callable from C.
In order to use the functions provided by a library, the user must provide the -F option on the command-line, followed by the name of the library (excluding filename extensions). An example library is included with the Vartools distribution, which may be used as a template.

Options
In addition to the commands mentioned above, a number of command-line options are available to change the behavior of Vartools. These options may be given at any point on the command-line (except for within the list of parameters supplied to a command). The available options include: -basename: only the base filename of each light curve is included in the output table, rather than the full path.
-bufferlines: use this to adjust the number of lines buffered within Vartools before being written to standard out during parallel processing. This buffering is separate from the system buffering done on stdout itself, which may be turned off with the -nobuffer option. Thebufferlines option only has an effect if the -parallel option is used, and will help speed up the processing by preventing threads from waiting on each other to output results. The tradeoff with using a larger buffer is that more memory is used, and results will be output less frequently.
-example: used in conjunction with the name of a command to print out an example of how that command is used.
-functionlist: shows the list of function names that are recognized in analytic expressions.
-header: include column headers in the output ascii table.
-headeronly: prints the output header that a call to Vartools would produce and exits without processing any light curves. This is useful for finding the column names or numbers to use when values computed by one command are used to set parameter values for subsequent commands.
-help: shows detailed help for a command.
-inlistvars: used to specify a format for the input light curve list, and store columns from this list as vector variables.
-inputlcformat: used to specify the format of the input light curves. One may use this option to read-in in more than just the default JD, magnitude, and magnitude uncertainty values from the light curve.
-jdtol: sets the tolerance for considering two observations to have come from the same epoch. The default tolerance is 10 −5 days. This is used by commands such as -TFA or -SYSREM which match points from different light curves.
-quiet: process the light curves, but do not output the ascii table of statistics.
-L: load a user-compiled library defining a new Vartools command.
-listcommands: shows a terse list of the available commands.
-log-command-line: includes the command-line in the header of the output ascii table.
-matchstringid: use image string identifiers, rather than the time, to match points from different light curves.
-nobuffer: outputs each line in the ascii table as soon as it is generated, rather than buffering the output.
-noskipempty: by default empty light curves are skipped and not included in the output table. To not skip these give this option. This option has no effect if the -readall option is used.
-numbercolumns: pre-pend column numbers to the header names in the output ascii table.
-oneline: display each computed statistic on a separate line for each light curve, rather than in a table format.
-parallel: process multiple light curves in parallel.
-randseed: specifies a seed for the random number generator. If this is not used, then every call to Vartools will produce the same set of random numbers. The system time may be used as a seed to produce quasi-random numbers.
-readall: read-in all of the light curves at once, rather than as they are processed.
-readformat: This option is deprecated in the current version of Vartools, users are encouraged to use theinputlcformat option instead. It is provided to allow compatibility for older processing scripts.
-redirectstats: output the ascii table of statistics to a file rather than to stdout.
-showinputlcformat: displays the expected format of the input light curves, given a call to Vartools. No processing will occur if this option is given, it is provided for assistance in developing scripts using Vartools.
-showinputlistformat: displays the expected format of the input list of light curves, given a call to Vartools. No processing will occur if this option is given, it is provided for assistance in developing scripts using Vartools.
-skipmissing: do not abort if a missing or unreadable light curve file is encountered. Instead skip the light curve and proceed with others in the list.
-tab: outputs the ascii table as a starbase table (tabdelimited columns).

Searching a set of light curves for periodic variable stars with Lomb-Scargle
For the first example, suppose we have a collection of light curves that we wish to search for periodicity using the L-S method. Further, suppose each light curve is stored in a separate file with the format shown in listing 4, and that we created a list of light curves as shown in listing 5. These light curves may be searched for periodic signals using the Generalized Lomb-Scargle periodogram as shown in listing 6. Here we search for periods between 0.01 days and 100.0 days sampling the periodogram for each light curve using frequency steps of 0.1/T , where T is the total time spanned by the light curve. We only output one peak, which is the highest one found, and we do not store the periodograms. The options "-header" and "-numbercolumns" are used to provide the header given in the output, and to also prepend the column number to the name of each column in the header.
As seen in the output the first light curve has a periodic signal of 77.7776 days, detected with false alarm probability 10 −5710 and spectroscopic S/N 387. In fact, this particular example light curve has a strong linear trend which is well-fit by a periodic signal with a period longer than the time-span of the data. The second light curve also has a periodic signal, this time with a period of 1.2344 days, detected with false alarm probability 10 −3999.6 . This example light curve has a strong sinusoidal signal ( Figure 1).
Listing 7 shows a command which will also output the periodograms for each light curve to the directory "PER" and identify 5 peaks in the periodogram, in each case whitening the light curve (i.e., fitting and removing a harmonic function from the light curve) at the previously identified peak before finding the next one. For clarity we truncate the output shown below, the actual output contains 21 columns, with the period, false alarm probability, periodogram value, and S/N reported for each of the five peaks. Figure 1 shows an example of the periodogram output for the second light curve, while the format of the periodogram file produced by Vartools is shown in listing 8.

Processing a FITS-format light curve from Kepler
Listing 9 shows an example of processing the quarter 1 binary FITS-format light curve from the NASA Kepler mission (Borucki et al., 2010) for the star KIC 1429092 11 . The -inputlcformat option indicates the table columns 11 This light curve may be downloaded from https://archive.stsci.edu/kepler/publiclightcurves.html to use. The first column is BJD − 2454833, the eigth column is the PDC-corrected simple aperture photometry flux, and the ninth column is the uncertainty on this flux. We use the -changevariable commands to tell Vartools to use the pdcsap f lux variable for the light curve brightness measurements, and the pdcsap f lux err variable for the uncertainties. Of course in this example we could have read columns 8 and 9 directly into the mag and err variables, but it can be helpful to keep track of what quantities are actually being input when there are different possibilities. The -fluxtomag command converts the flux into magnitudes. We then calculate the r.m.s. scatter of the magnitudes and the L-S periodogram. The magnitudeconverted light curve is output to the file "tmp.lc" using the specified format. Lines 32 through 35 of listing 9 show the format of the output light curve.

Running the BLS algorithm on a light curve to find
transits and eclipses Listing 10 shows an example of running the BLS algorithm on a light curve containing a transit signal. This example light curve is included with the Vartools distribution, and is produced in a manner similar to that described below in example 3.4.
The parameters are as follows: we consider fractional transit durations between 0.002 and 0.05 in phase, search for periods between 0.5 days and 100 days, use 100,000 frequency steps and 500 phase bins. We use 0 for the timezone (this is used to assign points to an individual night in calculating the fraction of ∆χ 2 that comes from one night), and report the top five peaks in the periodogram. We output the periodograms to the directory PER, output the best-fit transit model to the directory MODEL, and we do not subtract the model from the light curve. A trapezoid-shaped transit is fit to each peak to refine the transit parameters compared to those found in the full frequency scan. We use the "nobinnedrms" keyword to useSR and σ SR in equation 14 for the output BLS SN columns, rather thanSR(f) and σS R (see the discussion in section 2.2.1.6), and we output a uniformly sampled phase curve for the best-fit transit model to the directory MODEL, to use for plotting.
The results are shown in figure 2. The transit is recovered with a period of 2.1232 days, with a spectroscopic S/N of 4.4, and with S/N pink = 14.5 (Eq. 15). We note that the latter is probably the most useful of the metrics calculated here for selecting transit candidates.

Simulating and Recovering a Transiting Planet Light
Curve In this example we simulate a light curve with rednoise, add a transit signal to it, recover the transit with -BLS, and then fit a Mandel and Agol (2002) transit model to the light curve. A similar procedure might be repeated millions of times to determine the detection efficiency of a transit survey. The commands are shown in listing 11. Listing 8: Format of a periodogram file produced by the command in listing 7. We truncate the initial header (columns 6 through 11 repeat P(omega) and the logarithm of the FAP for whitening cycles 2 through 4).  Listing 9: Processing a light curve from the Kepler mission, as discussed in example 3.2 1 prompt > vartools -i k p l r 0 0 1 4 2 9 0 9 2 -2 0 0 9 1 6 6 0 4 3 2 5 7 _ l l c . f i t s \ -inputlcfor ma t t :1 ,pdcsap_flux :8 ,pdcsap_flu x_ e rr :9 \ 3 -changevaria bl e mag pdcsap_flux \ -changevar ia bl e err pdcsap_flux _e rr \   command, together with the model phase curve (output using the "ophcurve" keyword). Middle: The BLS spectrum using the spectroscopic S/N when the "nobinnedrms" keyword is given (i.e., this is equivalent to SR, with the average value subtracted, and divided by the standard deviation). Bottom: The BLS spectrum using the spectroscopic S/N when the "nobinnedrms" keyword is not given (Eq. 14). The primary difference is a change in scale due the difference betweenSR (Eq. 10) andSR (Eq. 13).
In this case only the times and magnitude uncertainties are read-in from the light curve file, which we use to define the time-base for our simulated observations, and we then initialize all magnitudes to a value of 10 using the -expr command.
We use the -addnoise command to add time-correlated noise to the light curve, adopting the exponential noise model (Eq. 70) with ρ = 0.05 days, σ r = 0.0015 mag, and σ w = 0.0015 mag. The -changerror command is called to set the magnitude uncertainties equal to the r.m.s. scatter of the light curve, and we use the -o command to output the simulated light curve to the file "MODEL/1.simnoise".
Next comes the -Injecttransit command which is used to add the transit signal to the light curve. Here we use the "fix" version of each parameter to specify the value to use on the command line, if one were to carry out a large transit injection/recovery simulation with Vartools, then the "list" or "rand" versions may be used. The injected signal has a period of 2.123456 days, a planet radius of 1.0 R J , a planet mass of 1.0 M J , a random phase, a random inclination, e and ω fixed to zero, a stellar mass of 1.0 M ⊙ , a stellar radius of 1.0 R ⊙ , and quadratic limb darkening coefficients of 0.3 and 0.2. We do not output the model, but use the following -o command to output the light curve with the transit injected to the file "MODEL/1.injecttransit".
The -BLS command is then used in a similar fashion to the previous example to recover the transit, though in this case we only identify 1 peak in the spectrum. We then use the -MandelAgolTransit command to fit a Mandel and Agol (2002) transit model to the light curve, initializing the parameters based on the previous -BLS command. We use the same limb darkening coefficients as for the injection, and vary all parameters except for e, ω, and the limb darkening coefficients. We do not attempt to fit an RV curve, and we do not subtract the best-fit model. We output the model together with the model phase curve to the directory "MODEL". The phase curve (for plotting) runs from phase −0.5 to 0.5 with a phase step of 0.001.
A selection of the output from this routine is shown in the listing. Figure 3 shows the simulated noise-only light curve, together with the phase-folded light curve with the transit signal and Mandel and Agol (2002) model overplotted.

Searching a set of light curves for variable stars, while
filtering noise In this example we perform a general search for variability on a set of light curves after applying two different noise-filtering techniques. Listing 12 shows the format of the light curves for this example. In addition to the usual time, magnitude and uncertainty values provided in the light curve, auxiliary information is provided for each observation including the filename of the image from which the observation is generated, the X and Y coordinates of the source on the image, the full-width at half maximum (FWHM), ellipticity and position angle of a Gaussian fit Listing 11: Simulating a light curve with red-noise, adding a transit signal, recovering it with BLS and fitting a Mandel and Agol (2002) transit model to it, as discussed in example 3.4 prompt > vartools -i EXAMPLES /1 \  B L S _ N p o i n t s i n t r a n s i t _ 1 _ 6 = 120 36 BLS_Ntransi t s_ 1_ 6 = 3 B L S _ N p o i n t s b e f o r e t r a n s i t _ 1 _ 6 = 169 38 B L S _ N p o i n t s a f t e r t r a n s i t _ 1 _ 6 = 69 BLS_Rednoi s e_ 1_ 6 = 0 .00113 40 BLS_Whiteno i se _ 1_ 6 = 0 .00221 B L S _ S i g n a l t o P i n k n o i s e _ 1 _ 6 = 17 .60859 42 B L S _ P e r i o d _ i n v t r a n s i t _ 6 = 6 .45330708 B L S _ d e l t a C h i 2 _ i n v t r a n s i t _ 6 = -327.05444 44 BLS_MeanMag_6 = 10 .00068 M a n d e l A g o l T r a n s i t _ P e r i o d _ 7 = 2 .12519667 46 M a n d e l A g o l T r a n s i t _ T 0 _ 7 = 53726 .95486357 MandelAg ol T r an s it _ r _7 = 0 .09941 48 MandelAgol T r an s it _ a _7 = 6 .80453 M a n d e l A g o l T r a n s i t _ b i m p a c t _ 7 = 0 .26979 50 M a n d e l A g o l T r a n s i t _ i n c _ 7 = 87 .72772 ... 52 M a n d e l A g o l T r a n s i t _ c h i 2 _ 7 = 0 .99461  to the point spread function (PSF), and the airmass of the observation.
Listing 13 shows the procedure. We will calculate five variability statistics: (1) the r.m.s. scatter of the light curve; (2) the Generalized L-S periodogram; (3) the BLS spectrum; (4) Stetson's J statistic; and (5) the harmonic AoV periodogram. This will be done three times for each light curve: on the raw (un-filtered) input light curve; on the light curve after decorrelating against the external parameters X, Y, FWHM, Ellipticity, PA and Airmass; and on the light curve after applying TFA filtering in addition to the decorrelation. To speed up the processing, this is run in parallel on the light curves.
Parameters are chosen for the commands as follows. For L-S we scan periods from 0.1 to 100 days with a subsample factor of 0.1. We only report the top peak, and we do not output the periodogram. For BLS we consider fractional transit durations between 0.002 and 0.05 in phase, search for periods between 0.5 days and 100 days, use 100,000 frequency steps and 500 phase bins. We use 0 for the timezone (this is used to assign points to an individual night in calculating the fraction of ∆χ 2 that comes from one night), and report the top five peaks in the periodogram. We do not output the periodogram or the best-fit transit model, and we do not subtract the model from the light curve. A trapezoid-shaped transit is fit to each peak to refine the transit parameters compared to those found in the full frequency scan. For Stetson's J we use a timescale of 0.5 days to distinguish between "near" and "far" observations, and "dates.txt" is a file listing all of the dates in the light curves. For the harmonic AoV we use a single harmonic model (i.e., a simple sinusoid), search for periods between 0.1 days and 100 days with a sub-sample factor of 0.1. A higher resolution scan, with a subsample of 0.01 is performed around each peak. We only report the top peak in the periodogram, and we do not output the periodogram. Note that the parameters chosen for this example are just to illustrate the use of the program, and are not meant to be best or default choices. In practice all parameters need to be carefully considered based on the properties of the light curves being analyzed.
After calculating the statistics for the raw light curve, we fit a model that is a linear function of all of the external parameters (in practice this is unlikely to be a very good model, one should inspect the light curves to determine the best function to use for filtering). This is then subtracted from the light curve using the 'correctlc' keyword. The linear coefficients which we optimize are labeled a0 Listing 13: Vartools procedure for carrying out a general search for variability on a set of light curves, including different levels of filtering, as discussed in example 3.5 1 prompt > vartools -l inputlist.txt \ -inputlcform at id :1 ,t :2 ,mag :3 ,err :4 ,x :5 ,y :6 ,fw :7 ,el :8 ,pa :9 ,air :10 \ ' a0 + a1 * x + a2 * y + a3 * fw + a4 * el + a5 * pa + a6 * air + a7 * sin (2 .0 * pi * x )+ a8 * sin (2 .0 * pi * y ) ' \ a 0 , a 1 , a 2 , a 3 , a 4 , a 5 , a 6 , a 7 , a 8 correctlc \ give a term that is periodic in one pixel. Light curves generated through aperture photometry commonly have systematic variations that are correlated with terms such as these. Following the '-linfit' call, we use '-expr' to set the light curve mean back to its pre-filtered value (subtracting the linear model causes the mean to be close to zero). This is achieved by adding the difference between the mean magnitude computed by '-rms' before '-linfit' and the mean magnitude computed after it. These are referenced by their output column names ("Mean Mag 0" and "Mean Mag 6"). The two calls to "-o" are used to output the light curves after each filtering step. The decorrelated light curves are output to the directory LCEPD, while the EPD+TFA-filtered light curves are output to the directory LCTFA.
To run TFA on the light curves we provide a list of files to be used as TFA templates in 'template.txt'. These light curves are assumed to have the same format as the input light curves, so we use the readformat option to indicate that column 1 will be used for matching points (this is the image id) and column 3 for the magnitudes. To use the image ids for matching points in the light curves to points in the templates, rather then the JD values, the '-matchstringid' option is given at the end of the call to Vartools. A list of all image ids is given in the file 'imagelist.txt'. Any template star that is within 10 pixels of the light curve being filtered will be removed from the template for that light curve. The xycol option is used to indicate that the average X and Y positions of the sources are read from columns 2 and 3 of the input light curve list file. The light curves passed on from this command will have the filtering applied, but we do not output the TFA coefficients or the model trend signals.
Running this command on a set of light curves produces an ascii table with 381 columns, which may be daunting to read. For convenience, Listing 14 shows a method for writing out the column names only, in a human-readable list. In this case one would fill in the '...' with the full command given in listing 13, and the '-headeronly' option causes Vartools to output the header and quit without processing any light curves. The calls to the standard command-line utilities 'sed' and 'awk' are used to reformat the names into a readable list, including blank lines between the columns produced by each command. We show the first few column names below the command. A simple selection of variables might be to select any light curve with log 10 FAP LS < −100, a BLS signal-to-pinknoise greater than 10, J stet > 1.5, or log 10 FAP AoV < −100. This can be done, for example, by applying a simple 'awk' script to the output. In practice these thresholds would need to be adjusted based on the noise properties of the survey (for example, the threshold of log 10 FAP LS < −100 is an extremely low false alarm probability which is a reasonable threshold for highly time-correlated noise, but far too conservative a threshold if the light curves actually Listing 14: A technique for printing out the column names for a Vartools procedure in a human-readable list. In this case the call to Vartools is the same as in listing 13, but with with the "-headeronly" option used to output only the command headers. The procedure is discussed in example 3.5. prompt > vartools -l inputlist.txt \ 2 -inputlcform at id :1 ,t :2 ,mag :3 ,err :4 ,x :5 ,y :6 ,fw :7 ,el :8 ,pa :9 ,air :10 \ -rms \ have Gaussian white noise).

Trend-filtering while preserving stellar variability with
signal-reconstruction TFA In this example we show the use of signal-reconstructionmode TFA (-TFA SR) to filter systematic noise from a light curve while preserving real variability. Here the variability signal is modeled using a harmonic series. The command is shown in listing 15. A periodic signal is first identified with the -LS command. We then save the light curve using the -savelc command, to allow us to try different filtering techniques in the same call to Vartools. The -Killharm command is used to subtract the signal from the light curve, and the call to -rms calculates the scatter of the residuals. Calling -restorelc 1 restores the light curve to its state at the -savelc command. We next run -TFA to apply the TFA method in non-reconstructive mode. The trends are read from the file "EXAMPLES/trendlist tfa", the list of dates from the file "EXAMPLES/dates tfa", we use a separation of 25.0 pixels to remove close neighbors from the trend list, we specify the columns in the input list from which to read the X and Y pixel positions of the source, and we subtract the model from the light curve but do not output the trend coefficients or model signal. After -TFA we output the filtered light curve, and then subtract the signal with -Killharm and calculate the scatter of the residuals with -rms. Finally we restore the light curve again to its original state, and then filter using -TFA SR. We use mostly the same options as for -TFA, but now we output the coefficients and model trends to the directory "EXAMPLES/OUTDIR1". We do not apply TFA first, we use a stopping threshold of 0.001 and a maximum number of iterations of 100 (in this case, where we are using a harmonic series to model the signal, the trend and model are fit simultaneously to the data so these parameters have no affect), we use a "harm"onic series to model the signal with zero higher-order harmonics and no sub-harmonics, and we take the period from the prior -LS command. Again we output the light curve, subtract the signal, and compute the scatter of the residuals.
By comparing the values for "RMS 4", "RMS 8" and "RMS 13" one can see the affect on the scatter in the residuals from the different filtering techniques. Applying -TFA in non-reconstructive model increases the scatter from 2.3 mmag to 7.2 mmag, but when the filtering is run in reconstructive mode the residual scatter decreases to 2.1 mmag. Figure 4 compares the phase-folded light curves after -TFA and after -TFA SR to the original light curve, where it can be seen that the signal-reconstruction correctly preserves the signal shape. Please note that this example is concocted primarily to illustrate the syntax and operation of the -TFA SR command. The noise-filtering in this example is not terribly effective, but in practice one would likely use a much larger set of template light curves, chosen in an intelligent fashion, which, if the light curves share systematic trends, would reduce the scatter of the residuals even further.
Listing 15: Example of filtering a light curve with TFA in non-reconstructive and reconstructive modes, as discussed in example 3.6 prompt > vartools -l EXAMPLES / lc_list_tf a _s r _h a rm -oneline -rms \ : Top: phase-folded light curve of a simulated variable star before applying noise filtering as discussed in example 3.6. This is the input light curve processed in listing 15. Middle: the phasefolded light curve after processing with the non-reconstructive-mode -TFA command (i.e., that output on line 8 of listing 15). Bottom: the phase-folded light curve after processing with the reconstructivemode -TFA SR command (i.e., that output on line 15 of listing 15). Using the -TFA SR command does a better job of filtering the noise while preserving the signal shape than does the -TFA command.

Using the Weighted Wavelet Z-Transform to Characterize a Quasi-Periodic Signal
Listing 16 shows an example of using the -wwz command to calculate the Weighted Wavelet Z-Transform of a light curve. The transform is calculated up to a maximum frequency of 2.0 d −1 with a frequency step of 0.25/T where T is the time base-line of the light curve. We use the "auto" keyword for both "tau0" and "tau1" to consider time-shifts running from the first observation to the last, and use a time-shift step of 0.1 days. The full transform is output to the file "EXAMPLES/OUTDIR1/8.wwz" while the projection of the wavelet onto the time-shift axis is output to the file "EXAMPLES/OUTDIR1/8.mwwz". These are plotted in Figure 5. We also show the light curve in Figure 6. For the full 2-d transform we use the "pm3d" keyword to output the file in a format that is convenient for plotting with the gnuplot program. The commands for plotting this are shown at the bottom of Listing 16. The data has a signal with a frequency of 0.3065 day −1 (P = 3.2632 days) centered at the time JD − 53725 = 10.174. The signal is not present at later times, while a lower significance signal with a frequency of ∼ 0.2 day −1 may be present at earlier times. This particular light curve comes from the MMT photometric survey of M37 presented by Hartman et al. (2008), and may be a spotted star with an evolving surface brightness distribution.

Performing an MCMC fit of a non-linear function to
a light curve Listing 17 shows an example of fitting an analytic function that is non-linear in its parameters to a light curve, and using an MCMC procedure to determine the uncertainties. The initial -stats command, together with the following -expr commands, are used to add a Gaussian function to a light curve. The -stats command is used to find the minimum and maximum times in the light curve, which are stored in the variables ST AT S t M IN 0 and ST AT S t M AX 0. The first two calls to -expr define the variables t1 and Dt in terms of the minimum and maximum values, and then the third call to -expr adds the Gaussian with amplitude 0.1 mag, standard deviation 0.05Dt and peak time t1+Dt * 0.2. For our illustrative call to -nonlinfit we provide the analytic function to fit to the light curve in terms of the free parameters a, b, c and d. We then set the initial values and uncertainty steps for each of these parameters on line 7 of listing 17. We use the "mcmc" keyword to indicate that the DEMCMC procedure should be run, stopping it after a total of 10000 links (in practice one would typically run more links, a relatively small number is used here to minimize the computing time for this simple example). The full MCMC chain will be output to the file "EXAMPLES/OUTDIR1/3.mcmc". The output table from Vartools gives the best-fit (maximum likelihood) values for the free parameters, together with the value of χ 2 for those parameters. Following this, the median and standard deviation of the MCMC chain for each parameter is provided. These can be taken as measures of the most likely values for each parameter, together with their 1σ uncertainties. Following the output from Vartools, lines 28 through 33 of listing 17 show the format of the file containing the full MCMC chain.
3.9. Converting the times in a light curve from UTC to BJD Listing 18 shows an example of converting the times in a light curve from UTC to TDB-corrected BJD. Lines 1-4 show the format of the input light curve file "EX-AMPLES/1.UTC", where the first column gives the UTC times of the observations. By default Vartools will attempt to read the time column directly into a doubleprecision variable, which it interprets as the JD. To convert the input time from UTC to JD, we use the call to -inputlcformat on line 7 of the listing. Here we indicate that the time will be read from the first column, that it is in "UTC" string format, and we then indicate how the string should be parsed, where "%Y" is interpretted as the year, "%M" is the month, "%D" is the day, "%h" is the UT hour, "%m" is the minute, and "%s" are the seconds. Finally we read the magnitude and uncertainties from columns 2 and 3, respectively. By default these are read in as double-precision variables, so we do not need to specify a variable type or format for these quantities. Lines 8 through 14 of the listing execute the time conversion. We first need to specify the input time system, which we do on line 8 (it is JD, and we use the "inputsys-utc" keyword to indicate that the input times have not been corrected for leap-seconds). We then indicate the output time system on line 9. Here we will convert to BJD, and will subtract 2400000 from the result, and we will convert to the TDB system to correct the times for leap-seconds. To convert to BJD we need to tell Vartools the coordinates of the star, which we do on line 10. We also need to provide a set of files from JPL which provide the data for calculating positions in the solar system ("ephemfile"), the list of leap-seconds that have occurred since 1970 ("leapsecfile") and the detailed position of locations on the surface of the Earth with respect to the center of the Earth-Moon system as a function of time ("planetdatafile") 12 . Line 14 indicates the observatory where the observations were made (in this case, Fred Lawrence Whipple Observatory in Arizona). Finally on line 15 we output the light curve to the file "EXAMPLES/OUTDIR1/1.bjdtdb", whose format may be seen on lines 17-20 of the listing.

Performance Tests
Here we present performance tests to compare the time required for I/O between ASCII and Binary format data, and to also determine the execution time for various periodfinding commands. We focus here on the periodic signal detection routines as these are generally the time-limiting components of processing pipelines for large datasets. 12 The relevant files may be obtained from ftp://naif.jpl.nasa.gov/pub/naif/generic_kernels/ Listing 17: Running an MCMC procedure to fit a function that is non-linear in its parameters to a light curve, as discussed in example 3.8.    The following tests were carried out on a machine with 32 AMD Opteron 6140 64-bit CPUs, each with a maximum clock speed of 2.6 GHz. The machine runs Debian GNU/Linux version 7.8 with Linux kernel version 3.2.0-4-amd64. Vartools was compiled on this machine with gcc version 4.7.2 using O2 optimization, and was linked against all supported external libraries. Here we are primarily interested in the processing time for different routines, so to minimize I/O overhead light curves were readin from shared memory (/run/shm/ in the listings below; this is a tmpfs file-system on this machine, and we confirmed that there was sufficient free memory during our tests that the files placed here were stored in RAM rather than swap space) and all output was ignored (redirected to /dev/null; we do not use the -quiet option, so time is still spent generating the ASCII table). Except where stated otherwise, all processing was performed on the light curve stored in "EXAMPLES/2.fits", which is a binary FITS table containing 3313 observations, and which is included in the Vartools distribution.

Test 1: ASCII vs. Binary I/O
There are four processes which contribute to the time required to input and output a light curve: (1) reading data from disk into memory; (2) parsing the raw data chunks into the appropriate variables within Vartools; (3) transforming data from the vector variables back into the appropriate format for storage on disk; (4) writing the data to disk. If the numerical data are stored in ASCII format then steps (2) and (3), which involve parsing ASCII strings into binary data and vice versa, can be time consuming. If the data is stored in binary format then steps (2) and (3) will be significantly faster (or not applicable, depending on the data organization). Storing the data in binary format also requires less total space than storing in ASCII to the same precision, reducing the time required for steps (1) and (4) as well. The primary benefit to using ASCII format is that the data are easily human-readable with widely available software tools, may be edited with simple text editors, and may be processed using a wide range of tools which expect ASCII format data (e.g., standard shell tools such as awk or sed).
To help in evaluating whether the time and storage savings from using binary data are worth the possible tradeoffs with data transparency, we carry out tests to evaluate the I/O performance of Vartools for handling both types of data.
The time required to read and write data from/to disk is highly system-dependent, and will vary in time as well, depending, for example, on whether the data have been recently accessed and are stored in cache memory. Regardless of the system, the time difference between ASCII and binary format data for steps (1) and (4) should be proportional to the difference in file-size. Here we focus on testing steps (2) and (3) and attempt to minimize the time required for steps (1) and (4) by using shared mem-ory for the file storage. We first generate an ASCII light curve using the command: Where the gauss() function used in the awk script is a simple routine for generating Gaussian random numbers, and where we use "%.17g" to write the ASCII data to full double-precision.
To determine the time for ASCII input only we carry out the following: Here the expression "-if 0 -rms -fi" is included to allow Vartools to execute without performing any processing of the light curve (if no commands are given to Vartools then it will terminate without reading the light curve, with this expression the light curve will be read-in and parsed into the appropriate variables, but the conditional placed in front of "-rms" will prevent it from executing). This test indicates that it takes Vartools ∼ 2.3 s of user+system time on our test machine to parse a 3-column ASCII table with 1,000,000 rows and with each column containing 17 digits. Repeated execution shows a variance of less than 0.1 s in the execution time. Note that running on an empty light curve takes less than 0.05 s so the initialization of Vartools is a negligible fraction of this time.
Using the "-l" option to execute on multiple copies of this same light curve shows that the time scales linearly with the number of light curves. Reducing the light curve size by half also reduces the time by half. If each of the input columns has 8 digit-precision rather than 17 digitprecision, then the process takes 1.5 s. If we use the "inputlcformat t:0,mag:0,err:0" option, in which case the light curve is read-from disk and the number of lines and columns are determined (this is done whenever Vartools reads in a light curve and requires a scan over all ASCII characters in the file), but the ASCII data are not converted to binary, then the process takes 1.2 s to execute. Finally, if we do also include the -rms command (i.e., the "-if 0" and "-fi" expressions are removed) the difference in time is negligible (less than the variance between repeated executions of the command), indicating that calculating the r.m.s. scatter of ASCII-format light curves is an I/O limited process, dominated by the ASCII to binary conversion, if not by the disk access itself.
To determine the time for ASCII output we carry out the following: The first command is used to convert the ASCII light curve to a binary FITS table. The second is used to generate a list with 100 repeated entries of the same light curve (this is done because we found that the input time for one of these light curves is not significantly longer than the time required to start Vartools), and the third performs the time test. Based on this test we find that it takes ∼ 0.1 s on this machine to read-in a FITS light curve with 1,000,000 rows and 3 double-precision columns (roughly 2/3 of the time is user time, the other third is system time). If we use the "-inputlcformat t:0,mag:0,err:0" option, then only 20 ms is required, with negligible system time.
To determine the output time for FITS we carried out the following test: Where the "nameformat" option is used to output the light curve to a different name than its input file. Subtracting the input time determined above we find that it takes ∼ 0.2 s to output the FITS light curve with 1,000,000 rows and 3 double-precision columns, or roughly twice as long as input.
Altogether we estimate that neglecting additional differences in disk access time when the data are stored on a physical hard disk, the input time for an ASCII light curve is approxiately 20 times longer than for a FITStable, while the output time is approximately 10 times longer for ASCII than for FITS. These scalings assume that the ASCII data consists of 3 columns each stored to 17-digit precision. Changing the precision, and/or the number of columns in the light curve will adjust the scalings. In terms of storage, the ASCII light curve with 1,000,000 data points takes 56 MB whereas the FITS-table with 3 double-precision columns takes 23 MB.
Note that based on applying similar tests to the file "EXAMPLES/2.fits", which is used for the performance tests discussed below, we conclude that it takes Vartools less than 0.3 ms to load the light curve.

Test 2: -LS
We now proceed to test the performance of individual Vartools commands. For our first test we executed the -LS command on 1000 copies of "EXAMPLES/2.fits". We choose the period range and subsample factors such that the periodogram would be evaluated at 16,384 frequencies (due to the use of FFTs within the code the number of frequencies sampled internally increases in powers of 2, and thus the execuction time shows discrete jumps as ceil(log 2 (N f )) increments). Tests were carried out with commands similar to the following: Repeated execution shows a variance of a few seconds in the time. The command takes on average 330 ms of user+system time per light curve on this system. When executed in parallel on 10 processors the routine takes 38 ms per light curve processed. Note that here, and in the examples below, we are calculating the periodogram for 10 different light curves in parallel, parallel processing of a single light curve is not currently supported. The total time scales close to linearly with the number of processors, however shared overhead and locks used to prevent conflicts between threads causes the scaling to be not exactly linear. In general the scaling is closer to linear with the number of processors when the execution time for a single process is longer. When 2048, 4096 or 8192 frequencies are scanned the process takes 20 ms, 48 ms, and 110 ms per light curve, respectively. Adjusting the number of points in the light curve, while keeping the number of frequencies fixed, has a negligible effect on the processing time which is dominated by the calculation of FFTs over the sampled frequencies. For very large light curves, however, there may be a difference as the time required for the "extirpolation" routine, which depends on the number of points in the light curve, becomes non-negligible.

Test 3: -BLS
We performed similar tests on -BLS to those conducted for -LS, using commands such as: The relevant numbers controlling the execution time in this case are the number of frequencies sampled N f (the above listing uses 10,000), the number of phase bins N bin (200 above), the maximum number of phase bins for the transit duration N bin,tr (for "q 0.005 0.05" and N bin = 200 above we have N bin,tr = 0.05 × 200 = 10), and the number of observations N obs . We expect the execution time to scale roughly as for some constants a and b. For the above example we find that the execution takes 1.1 s per light curve, or 120 ms of user+system time per light curve if we use 10 parallel processes (if the "nobinnedrms" keyword is used, then the execution is 920 ms per light curve, or 93 ms when running with 10 parallel processes; the difference in time is due to an extra square root being included in the constant b when the keyword is not used). Using 100,000 frequencies increases the execution time by a factor of ten as expected (11 s per light curve), while doubling N bin increases the execution time by a factor of approximately two (or 1.5 for "nobinnedrms").

Test 4: -aov harm
We performed similar tests on -aov harm to those conducted for the other commands. This was done using commands such as: Here we are scanning the same number of frequencies as for -LS and using a pure sinusoid for the signal. The processing however is signficantly slower than for -LS, with 16,384 frequencies taking 9.3 s per light curve when not running in parallel (or 930 ms per light curve when running with 10 parallel processes). Using half as many frequencies reduces the processing time by a factor of two. Using 2 or 3 harmonics (i.e.,the fundamental and first higher-order harmonic in the first case, and the fundamental and first two higher-order harmonics in the second case) increases the execution time to 11 s per light curve and 13 s per light curve, respectively (when not running in parallel).

Test 5: -wwz
We used the following command to test the performance of -wwz: Here we are scanning 249 frequencies each at 311 timeshifts. With 10 parallel processes this takes 1.3 s per light curve, or 13 s per light curve per processor. Doubling the number of frequencies sampled increases the execution time by a factor of two, as does doubling the number of time-shifts.

Test 6: -fastchi2
We used the following command to test the performance of -fastchi2 which is included as an extension to Vartools: where we explicitly load the library with the -L option to avoid spending time performing a disk search for the appropriate library, and the parameters for -fastchi2 are chosen to calculate 16,384 frequencies using a simple sinusoid. We find that this process takes 54 ms per light curve, or 5.6 ms per light curve if run on 10 parallel processors. When 2048, 4096 or 8192 frequencies are scanned the process takes 34 ms, 36 ms and 42 ms, respectively. Thus for large number of frequencies the -fastchi2 procedure is substantially faster than -LS (330 ms for 16,384 frequencies), but for 2048 or fewer frequencies the -LS command is faster to execute (20 ms per light curve).

Future Development
There are several areas in which Vartools may be improved to make it more useful for astronomical time series analysis. Several new processing commands are being actively developed, such as Fourier filtering methods, and techniques for stitching together different light curves for a given source. Other topics for future development, which require more fundamental changes to the organization of the code, include the parallelization of individual commands (currently only parallel processing of multiple light curves is supported), and support for parallel operation in non-shared-memory mode (e.g., using MPICH). The methods for handling analytic expressions may be extended to support vector and matrix operations. We are also working on support for user-developed commands and functions written in python. As a separate effort we are also working on developing a Graphical User Interface (GUI) for Vartools, written in python, which may be used to construct processing pipelines and interact with the results.

Appendix A. Command Syntax
Here we describe in detail the expected syntax for each command, including a brief explanation of the input parameters. Commands are listed in alphabetical order. See section 2.2 for a description of the algorithms and uses of each command. In general a command is executed by typing the name of the command (including the "-" prefix, e.g., "-rms" for the command which calculates the r.m.s. scatter in a light curve) followed by various parameters and options. Terms within angular brackets "< >" are required parameters. Terms within square brackets "[ ]" are optional parameters that may be ignored. Within a given set of brackets the pipe "|" is used to distinguish between different choices that are allowed (in general only one of the possible choices may be used). Any term not placed within brackets is required. Terms within quotations indicate keywords that should be typed out exactly as they appear, but without the quotation marks (these are typically used to control optional behavior, or to specify the name of a parameter which then follows). The keywords are case sensitive. Unquoted terms represent parameters for which the user should substitute an appropriate value when calling the command. Vartools obeys a strict ordering of parameters, you must provide them in the exact order indicated. Parameters referred to below as "flags" are used to control the behavior of the command. In these cases the user either substitutes "1" for the flag to turn the option on, or "0" to turn it off.
Many of the commands share syntax for controlling how parameter values should be determined. The syntax 1 <" f i x " v a l | " l i s t " [ " column " c o l ] | " f i x c o l u m n " <colname | colnum>> is adopted to allow the user to choose between fixing a parameter to a specific value for all light curves (the user types "fix" and then substitutes the appropriate value for val), to read the parameter value from the input list of light curve files allowing a different value for each light curve (the user types "list" then optionally indicates which column the parameter is to be read-in from by typing "column" and the column number; if the "column" keyword is not used, then the values are read-in from the next unused column in the file), or to use a value computed by a previously executed command (the user types "fixcolumn" and then gives either the name of the column, which is the column heading when the "-header" option is used, or the number of the column in the output table). In some cases only a subset of these options is allowed. Below we list the expected syntax for each command, followed by a table explaining the various parameters.
1 −addnoise < " w h i te " 3 <" s i g w h i t e " <" f i x " v a l | " l i s t " [ " column " c o l ]>> | " s q u a r e e x p " 5 <"rho " <" f i x " v a l | " l i s t " [ " column " c o l ]>> <" s i g r e d " <" f i x " v a l | " l i s t " [ " column " c o l ]>> 7 <" s i g w h i t e " <" f i x " v a l | " l i s t " [ " column " c o l ]>> [ " b i n t i m e " <" f i x " v a l | " l i s t " [ " column " c o l ] > ] 9 | " exp " <"rho " <" f i x " v a l | " l i s t " [ " column " c o l ]>> 11 <" s i g r e d " <" f i x " v a l | " l i s t " [ " column " c o l ]>> <" s i g w h i t e " <" f i x " v a l | " l i s t " [ " column " c o l ]>>

13
[ " b i n t i m e " <" f i x " v a l | " l i s t " [ " column " c o l ] > ] | " matern " 15 <"nu" <" f i x " v a l | " l i s t " [ " column " c o l ]>> <"rho " <" f i x " v a l | " l i s t " [ " column " c o l ]>> 17 <" s i g r e d " <" f i x " v a l | " l i s t " [ " column " c o l ]>> <" s i g w h i t e " <" f i x " v a l | " l i s t " [ " column " c o l ]>>

19
[ " b i n t i m e " <" f i x " v a l | " l i s t " [ " column " c o l ] > ] | " w a v e l e t " 21 <"gamma" <" f i x " v a l | " l i s t " [ " column " c o l ]>> <" s i g r e d " <" f i x " v a l | " l i s t " [ " column " c o l ]>> 23 <" s i g w h i t e " <" f i x " v a l | " l i s t " [ " column " c o l ]>> >   − The periodogram will be sampled at a resolution of subsample/T where T is the time baseline of the light curve. finetune − Peaks found in the initial scan will then be resampled at a resolution of finetune/T . Npeaks − The number of peaks in the periodogram to report. operiodogram − Flag indicating if the periodogram for each light curve should be output to a separate file. outdir − Required if operiodogram is set. The periodogram will be output to the file $outdir/$basename.aov where $outdir is the value specified for outdir, and $basename is the base filename of the light curve (any leading directory is stripped). "whiten" − Keyword used to whiten the light curve and recalculate the periodogram for each peak. By default this is not done. "clip" − Adjust the clipping performed on the periodogram in calculating the spectroscopic S/N of a peak. By default iterative 5σ clipping is performed. clip − The sigma-clipping factor to use. clipiter − Flag which if set causes iterative clipping to be performed. "uselog" − Keyword used to calculate the S/N from the natural logarithm of the AoV spectrum. By default the S/N is calculated directly from AoV. "fixperiodSNR" − Report also AoV, the false alarm probability and the S/N for a fixed period. "aov" − Use the peak period from the last executed -aov command for the fixperiodSNR. "ls" − Use the peak period from the last executed -LS command. "injectharm" − Use the period from the last executed -Injectharm command (the first one is used if multiple periods in a given command. − The periodogram will be sampled at a resolution of subsample/T where T is the time baseline of the light curve. finetune − Peaks found in the initial scan will then be resampled at a resolution of finetune/T . Npeaks − The number of peaks in the periodogram to report. operiodogram − Flag indicating if the periodogram for each light curve should be output to a separate file. outdir − Required if operiodogram is set. The periodogram will be output to the file $outdir/$basename.aov harm where $outdir is the value specified for outdir, and $basename is the base filename of the light curve (any leading directory is stripped). "whiten" − Keyword used to whiten the light curve and recalculate the periodogram for each peak. By default this is not done. "clip" − Adjust the clipping performed on the periodogram in calculating the spectroscopic S/N of a peak. By default iterative 5σ clipping is performed. clip − The sigma-clipping factor to use. clipiter − Flag which if set causes iterative clipping to be performed. "fixperiodSNR" − Report also AoV, the false alarm probability and the S/N for a fixed period. "aov" − Use the peak period from the last executed -aov command for the fixperiodSNR. "ls" − Use the peak period from the last executed -LS command. "injectharm" − Use the period from the last executed -Injectharm command (the first one is used if multiple periods in a given command. The files will be written to $outdir/$basename.autocorr where $outdir is the value specified for outdir, and $basename is the base filename of the light curve (any leading directory is stripped).

− a u t o c o r r e l a t i o n s t a r t s t o p s t e p o u t d i r
1 − b i n l c <"a v e r a g e " | " median " | " w e i g h t e d a v e r a g e "> <" b i n s i z e " b i n s i z e | " n b i n s " n b i n s> <" t c e n t e r " | " t a v e r a g e " | " tmedian " | " t n o s h r i n k " [ " b i n c o l u m n s o n l y "] > Table A.5: Input Parameters for "-binlc" Command (Section 2.2.6.1) "average" − Take the average of points in a bin "median" − Take the median of points in a bin "weightedaverage" − Take the weighted average of points in a bin "binsize" − Specify the size of the bin in the time unit of the light curve. Type the "binsize" keyword and then provide the value. "nbins" − Specify the number of bins to use. In this case the binsize is equal to T /nbins where T is the time-span of the light curve. Type the "nbins" keyword and then provide the value. "bincolumns" − By default all light curve vectors except for the time, and uncertainty, will be binned using the "average", "median", or "weightedaverage" as specified by the user. Any vectors that the user wishes to bin in a different manner may be specified with this option. Following the "bincolumns" keyword provide a comma separated list of light curve variables and associated statistics. After a given variable use a ":" and then indicate the statistic to calculate. The options for the statistics are the same as for the -stats command.
"firstbinshift" − By default the first bin begins at the initial time in the light curve (t 0 ). Use this keyword to change this behavior, in which case the first bin will start at t 0 − firstbinshift/binsize. "tcenter" − Take the binned light curve time to be the time at the center of each bin. "taverage" − Take the binned light curve time to be the average of the times of points that fall in a given bin. "tmedian" − Take the binned light curve time to be the median of the times of points that fall in a given bin. "tnoshrink" − Do not shrink the light curve. In this case all points in the light curve will be replaced by their binned value, but the times will not be changed. "bincolumnsonly" − If this keyword is given, then only the columns specified after the "bincolumns" keyword will be binned. All other light curve vectors will be left unchanged. This option is only available when the "tnoshrink" option is used.  − Allow for a period dependent limit on the transit durations to search. r corresponds approximately to the stellar radius for lower main sequence stars assuming central transits, circular orbits, and that times are specified in days. rmin − The minimum transit duration at period P will be 0.076rmin 2/3 P 1/3 . rmax − The maximum transit duration at period P will be 0.076rmax 2/3 P 1/3 . "q" − Use a period-independent limit on the fractional transit durations considered. qmin − The minimum transit duration at period P will be qminP . qmax − The maximum transit duration at period P will be qmaxP . minper − The minimum period to search for transits. maxper − The maximum period to search for transits. nfreq − The number of frequencies to search (the frequency step-size is ∆f = (1/minper − 1/maxper)/nfreq). nbins − The number of phase bins to use. timezone − The timezone of the observatory (in hours from UTC). This only affects the grouping of points used to determine the fraction of ∆χ 2 from a single day. If multiple observatories were used, or if the observations were performed from space one can simply provide 0 for this value and ignore the reported ∆χ 2 from one day. Npeak − The number of peaks to identify in the BLS spectrum. outperiodogram − A flag used to indicate whether or not the BLS spectra should be output. outdir − Required if outperiodogram is set to 1. The BLS spectra will be written to files with the name $outdir/$basename.bls $outdir is the value specified for outdir, and $basename is the base filename of the light curve (any leading directory is stripped). omodel − A flag used to indicate whether or not the best-fit transit model light curves should be output. These light curves are evaluated at the observed times. model outdir − Required if omodel is set to 1. The models files will be written to $model outdir/$basename.bls.model. correctlc − A flag used to indicate whether or not the best-fit transit model should be subtracted from the observations before passing the light curve on to the next command. "fittrap" − An optional keyword which, if used, causes a trapezoid-shape transit model to be fit to the light curve around each BLS peak identified. "nobinnedrms" − An optional keyword used to adjust the way in which the BLS spectroscopic S/N statistic is calculated. If this keyword is used then the average and standard deviation of the spectrum are calculated only from the maximum SR values at each trial frequency. If it is not used, then these are calculated from all SR values considered (including transit epochs and durations that do not optimize the fit at a given frequency). Using this keyword speeds up the BLS calculation, but the S/N will tend to be suppressed for high significance detections.
"ophcurve" − A keyword used to output best-fit model light curves that are uniformly sampled in phase. outdir − The models will be written out to files named $outdir/$basename.bls.phcurve. phmin − The starting phase to use in the model. phmax − The ending phase to use in the model. phstep − The phase step-size to use in the model. "ojdcurve" − A keyword used to output best-fit model light curves that are uniformly sampled in time between the first and last observed times in the light curve. outdir − The models will be written out to files named $outdir/$basename.bls.jdcurve. jdstep − The time step-size to use in the model. "stepP" − Sample the BLS spectrum at uniform steps in period (by default uniform steps in frequency are used). "steplogP" − Sample the BLS spectrum at uniform steps in the logarithm of the period (by default uniform steps in frequency are used). "adjust-qmin-bymindt" − Adaptively set the minimum q value to the maximum of qmin or ∆t min × f where ∆t min is the minimum time difference between consecutive points in the light curve. "reduce-bins" − Adaptively reduce the number of phase bins at each frequency such that there are no more than two bins to sample a transit of phase duration qmin. | " e x p r " expr> <"r " rmin rmax | " q" qmin qmax >  Table A.7: Input Parameters for "-BLSFixPer" Command (Section 2.2.1.8) "aov" | "ls" | "list" | "fix" | "fixcolumn" | "expr"

−
Source for the period at which to calculate the BLS model (either from the last -aov command, the last -LS command, from the input light curve list file, fixed to a value given on the command line, set to the output from a previously executed command, or determined by evaluating an analytic expression). "r" − Allow for a period dependent limit on the transit durations to search. r corresponds approximately to the stellar radius for lower main sequence stars assuming central transits, circular orbits, and that times are specified in days. rmin − The minimum transit duration at period P will be 0.076rmin 2/3 P 1/3 . rmax − The maximum transit duration at period P will be 0.076rmax 2/3 P 1/3 . "q" − Use a period-independent limit on the fractional transit durations considered. qmin − The minimum transit duration at period P will be qminP . qmax − The maximum transit duration at period P will be qmaxP . nbins − The number of phase bins to use. timezone − The timezone of the observatory (in hours from UTC). This only affects the grouping of points used to determine the fraction of ∆χ 2 from a single day. If multiple observatories were used, or if the observations were performed from space one can simply provide 0 for this value and ignore the reported ∆χ 2 from one day. omodel − A flag used to indicate whether or not the best-fit transit model light curves should be output. These light curves are evaluated at the observed times. model outdir − Required if omodel is set to 1. The models files will be written to $model outdir/$basename.blsfixper.model. correctlc − A flag used to indicate whether or not the best-fit transit model should be subtracted from the observations before passing the light curve on to the next command. "fittrap" − An optional keyword which, if used, causes a trapezoid-shape transit model to be fit to the light curve around each BLS peak identified.  − Optional keyword used to indicate that the transit depth should be fixed (by default it is allowed to vary). If this keyword is used, then one should follow it with an indication of how the fixed depth should be determined. "qgress" − Optional keyword used to indicate that the fractional duration of transit ingress (i.e., the ingress duration divided by the total transit duration) should be fixed (by default it is allowed to vary). This keyword may only be used if "fixdepth" is also used. One should follow the keyword with an indication of how the fixed fractional ingress duration should be determined. minper − The minimum period to search for transits. maxper − The maximum period to search for transits. nfreq − The number of frequencies to search (the frequency step-size is ∆f = (1/minper − 1/maxper)/nfreq). timezone − The timezone of the observatory (in hours from UTC). This only affects the grouping of points used to determine the fraction of ∆χ 2 from a single day. If multiple observatories were used, or if the observations were performed from space one can simply provide 0 for this value and ignore the reported ∆χ 2 from one day. Npeak − The number of peaks to identify in the BLS spectrum. outperiodogram − A flag used to indicate whether or not the BLS spectra should be output. outdir − Required if outperiodogram is set to 1. The BLS spectra will be written to files with the name $outdir/$basename.blsfixdurtc $outdir is the value specified for outdir, and $basename is the base filename of the light curve (any leading directory is stripped). omodel − A flag used to indicate whether or not the best-fit transit model light curves should be output. These light curves are evaluated at the observed times. model outdir − Required if omodel is set to 1.
The models files will be written to $model outdir/$basename.blsfixdurtc.model. correctlc − A flag used to indicate whether or not the best-fit transit model should be subtracted from the observations before passing the light curve on to the next command. "fittrap" − An optional keyword which, if used, causes a trapezoid-shape transit model to be fit to the light curve around each BLS peak identified. "ophcurve" − A keyword used to output best-fit model light curves that are uniformly sampled in phase. outdir − The models will be written out to files named $outdir/$basename.blsfixdurtc.phcurve. phmin − The starting phase to use in the model. phmax − The ending phase to use in the model. phstep − The phase step-size to use in the model. "ojdcurve" − A keyword used to output best-fit model light curves that are uniformly sampled in time between the first and last observed times in the light curve. outdir − The models will be written out to files named $outdir/$basename.blsfixdurtc.jdcurve. jdstep − The time step-size to use in the model.
− c h a n g e e r r o r No parameters or options are available.
1 − c h a n g e v a r i a b l e <"t " | "mag" | " e r r " | " i d"> v a r    The σ-clipping factor to use. If a value <= 0 is specified, then σ-clipping is not performed, but points with errors <= 0 or NaN magnitude values will be clipped from the light curve. iter − A flag that is 1 for iterative clipping (performed continuously until no further points are removed), or 0 to not do continuous iterative clipping. "niter" − An optional keyword used to specify a fixed number of iterations to use in performing the clipping. "median" − By default clipping is done with respect to the mean. Use this keyword to cause the clipping to be done instead with respect to the median.

− c l i p s i g c l i p i t e r [ " n i t e r " n ] [ " median " ]
1 − c o n v e r tti m e <"i n p u t " <"mjd" | " j d " | " h j d " | " b j d " >>

3
[ " i n p u t s u b t r a c t " v a l u e ] [ " i n p u ts y s − td b " | " i n p u t s y s − u t c " ] <"o u tp u t " <"mjd" | " j d " | " h j d " | " b j d " >>  Table A.12: Input Parameters for "-converttime" Command (Section 2.2.6.3) "input" − Indicate the input time system used. The required "input" keyword must be followed by one of the following keywords. "mjd" − Modified Julian Date (MJD = JD − 2400000.5). "jd" − Julian Date. "hjd" − Helio-centric Julian Date. "bjd" − Bary-centric Julian Date. "inputsubtract" − Optional keyword used to indicate a constant that has been subtracted from the input times (i.e., if the input time is HJD − 2400000, one would give "input hjd inputsubtract 2400000" on the command-line). "inputsys-tdb" − Optional keyword used to indicate that the input times have been corrected for leap-seconds. "inputsys-utc" − Optional keyword used to indicate that the input times have not been corrected for leap-seconds (this is assumed by default). "output" − Indicate the output time system used. The required "output" keyword must be followed by a keyword indicating the time system, options are the same as for "input". "outputsubtract" − Optional keyword used to indicate a constant that is to be subtracted from the output times. "outputsys-tdb" − Optional keyword used to indicate that the output times should be on a system corrected for leap-seconds. "outputsys-utc" − Optional keyword used to indicate that the output times should be on a system that is not corrected for leap-seconds (by default the same system as the input times is assumed). "radec" − Indicate the source for the RA and Dec coordinates of the source. This is required if converting to/from HJD or BJD. Both coordinates should be in decimal degrees. "input-radec" − Optionally specify a different soruce for the RA and Dec coordinates that were assumed in calculating HJD or BJD for the input times. This may be useful, for example, if the times were converted to HJD or BJD assuming the coordinates for the center of an image, but one wishes now to correct them to the coordinates of the source itself. "epoch" − Optionally specify a time epoch in years on the Common Era system for which the input RA and Dec coordinates are provided. The default is 2000.0. "input-ppm" − Optionally indicate the source for the input RA and Dec proper motions assumed for the source. These should be given in units of milliarcseconds per year. By the initial times are assumed to have been calculated with no proper motion, or with the proper motion specified using the "ppm" keyword. "ephemfile" − Provide the filename for the JPL NAIF ephemeris file used to determine the position of the Earth with respect to the Solar System Barycenter. If not specified Vartools will check if the environment variable CSPICE EPHEM FILE has been set, and use it if it has been set. "leapsecfile" − Provide the filename for the JPL NAIF leap-second file used to determine the number of leap-seconds since a reference epoch. This can also be set with the environment variable CSPICE LEAPSEC FILE. "planetdatafile" − Provide the filename for the JPL NAIF planetary physical data file used to determine the location of the observer with respect to the center of the Earth in the J2000.0 intertial frame. This can also be set with the environment variable CSPICE PLANETDATA FILE. "observatory" − Optionally provide this keyword to specify an observatory used to make the observations. If this is set then one should either give: code − A short string used to indicate the observatory, e.g., "maunakea" for Mauna Kea Observatory, Hawaii, "lasilla" for La Silla Observatory, etc. "show-codes" − If this keyword is given, then the list of allowed observatory codes will be output, and Vartools will not perform any processing. "coords" − Optionally provide the latitude (in degrees), longitude (in degrees east of the prime meridian), and altitude (in meters) of the observatory at which the measurements were performed. One must then indicate the source for these quantities. In addition to the standard "fix" and "list" keywords, one may also use the "fromlc" keyword to read these quantities as columns in the light curve. This may be used if observations carried out from different locations were combined into a single light curve file.
1 −copylc N c o p i e s   Table A.14: Input Parameters for "-decorr" Command (Section 2.2.3.1) correctlc − A flag indicating whether or not the light curves passed on to the next command will be decorrelated. If it is set to 0 then the value of χ 2 from the decorrelation together with the decorrelation coefficients will be output to the   The number of points per 1/T frequency element to include in the calculated power spectrum, where T is the time base-line of the input light curve. "maxfreq" − An optional keyword which should be followed by the maximum frequency (in cycles per day) at which to calculate the power spectrum. If not specified, the power spectrum will be calculated up to the Nyquist frequency. "outdspec" − An optional keyword used to output the raw Discrete Fourier Transform (i.e., before applying the CLEAN algorithm) to a file. dspec outdir − Required if the "outdspec" keyword is set. The raw DFT will be output to a file named $dspec outdir/$basename.dftclean.dspec where $basename is the base filename of the input light curve (stripped of any leading directory names). "finddirtypeaks" − An optional keyword used to search the raw DFT for peaks. Npeaks − The number of peaks to report from the raw DFT power spectrum. This parameter is required if the "finddirtypeaks" keyword is given. "clip" − Optional keyword to change the clipping method used in calculating the average and standard deviation of the raw power spectrum for determining the spectroscopic S/N of each peak. If this keyword is given then follow it with the σ-clipping factor and a flag to indicate whether or not continuous iterative clipping should be performed. By default iterative 5σ clipping is performed. "outwfunc" − An optional keyword used to output the window function to the file $wfunc outdir/$basename.dftclean.wfunc. "clean" − If this keyword is specified then the CLEAN deconvolution algorithm will be executed on the power spectrum. gain − The gain factor used in cleaning the power spectrum. It should have a value between 0.1 and 1. SNlimit − Cleaning will proceed until the last peak has a value that is less than SNlimit times the standard deviation. outcbeam − Optional keyword to output the "clean beam" used in producing the final power spectrum. This will be output to the file named $cbeam outdir/$basename.dftclean.cbeam. findcleanpeaks − An optional keyword used to search the cleaned DFT for peaks. It takes the same additional parameters as "finddirtypeaks". "useampspec" − Give this keyword to use the amplitude spectrum for calculating the S/N of peaks rather than the power spectrum, which is used by default. "verboseout" − Give this keyword to output the average and standard deviation of the spectrum before and after clipping, in additional to the final S/N value.  − An optional keyword used to specify the column from the input light curve list with the reference magnitude of the source. If not given, then the next available column in the light curve list file will be used.
1 − e n s e m b l e r e s c a l e s i g s i g c l i p    − If this keyword is specified then matchrad is assumed to be in arcseconds and the X and Y positions of the sources are RA and Dec. in decimal degrees. If this keyword is not given, then rectangular matching will be performed on the X and Y coordinates, and matchrad should be in the same units. "xycol" − An optional keyword for specifying the columns from the input light curve list file to use for determining the X and Y coordinates of the each of the light curves. If not given, then the next unused columns in the list file will be assumed. "fix" | "list" | "fixcolumn" − Source for the period to use for the light curve.
"starlist" − By default the input light curve list is matched to itself. Use this keyword to match instead to the list in the file starlistfile. The file should have 3 white-space delimited columns, the first being a string identifier, the second and third being the X and Y coordinates of the source. "zeromag" − Optional keyword for specifying the zero-point magnitude for converting from magnitudes into fluxes. The default value is 25.0. "nofluxconvert" − If this optional keyword is given then the magnitude to flux conversion is not performed. "Nharm" − Optional keyword to change the number of harmonics to include in the Fourier series fit to each light curve to determine its amplitude of variability. If it is 0 then only a sinusoid will be fit. The default value is 2. "omatches" − Optional keyword to output the names and flux amplitudes of all stars matching to each potential variable. outputmatchfile is the name of the file to output this information to.
−fluxtomag m a g c o n s ta n t o f f s e t Table A.20: Input Parameters for "-fluxtomag" Command (Section 2.2.6.7) mag constant − The magnitude of a source with a flux of 1 ADU. offset − A constant to add to the output light curve magnitudes. Table A.21: Input Parameters for "-GetLSAmpThresh" Command (Section 2.2.1.2) "ls" | "list" − The source for the period for which to calculate the minimally detectable amplitude. If the keyword "ls" is given then the source is the highest peak found by the most recent -LS command. If it is "list" then the period is read-in from the input light curve list. minp − The minimum period considered by the L-S run. This is needed for determining the false alarm probability, which depends on the period range scanned. thresh − The maximum value of log 10 FAP that the light curve have had and still have been considered a detection. "harm" − Use this keyword to fit a harmonic series to the light curve to determine its amplitude. Nharm − The number of harmonics to include in the series. Nsubharm − The number of sub-harmonics to include in the series. "file" − Use this keyword to read-in the model signal from a file (in general a separate file is used for each light curve).
listfile − The name of a file listing the files containing the model signals for each of the light curves. This list file has two columns of the form: signal file signal amp, with one line for each light curve being processed. Each signal file should contain the signal magnitude in the third column. signal amp is the amplitude in magnitudes of the signal (to allow for cases where the signal amplitude is greater than the difference between the minimum and maximum values in the file). "noGLS" − If this optional keyword is given then the traditional L-S periodogram is used for determining the false alarm probability of a signal. By default the Generalized L-S is used.  − An analytic expression evaluated for a given light curve. If this evaulates to a number different from 0 the expression will be treated as "true." -command1 ...
-commandN − A set of commands that will be executed conditional upon testing the expression.
−Injectharm <" l i s t " [ " column " c o l ] | " f i x " p e r 2 | " rand " minp maxp | " l o g r a n d " minp maxp | " r a n d f r e q " minf maxf 4 | " l o g r a n d f r e q " minf maxf> Nharm (<" a m p l i s t " [ " column " c o l ] 6 | " ampfix " amp | " amprand" minamp maxamp | " amplogrand " minamp maxamp> [ " amprel " ]  − Take the period of the harmonic signal to add to the light curve from the input light curve list file, optionally specifying the column number. "fix" − Fix the period to the specified value. "rand" − Choose a random number for the period between minp and maxp. "logrand" − Choose a random number for the period from a distribution that is uniform in logarithm between minp and maxp. "randfreq" − Choose a random number for the period from a distribution that is uniform in frequency between minf and maxf (in cycles per day). "lograndfreq" − Choose a random number for the period from a distribution that is uniform in the logarithm of the frequency between minf and maxf (in cycles per day). Nharm − The number of harmonics to include in the harmonic series. If this is 0 then only the fundamental mode is included (i.e., a simple sine curve). For each of the Nharm + 1 modes one must specify the source for the amplitude and phase. This done with the following keywords: "amplist" − Take the amplitude from the input light curve list file. Optionally specify the column to use with the "column" keyword. "ampfix" − Fix the amplitude to the specified value. "amprand" − Use a random number drawn from a uniform distribution between minamp and maxamp. "amplogrand" − Use a random number drawn from a uniform-log distribution between minamp and maxamp. "amprel" − If specified, amplitudes for harmonics 1 through Nharm are given relative to the amplitude of the fundamental mode.
"phaselist" − Take the phase from the input light curve list file. Optionally specify the column to use with the "column" keyword. "phasefix" − Fix the phase to the specified value. "phaserand" − Draw a random number between 0 and 1. "phaserel" − If specified, phases for harmonics 1 through Nharm are given relative to the fundamental mode. Nsubharm − The number of subharmonics to use. For each of the Nsubharm subharmonics one must specify the source for the amplitude and phase as above, the "amprel" and "phaserel" keywords may also be provided to make the amplitude and/or phase relative to the values for the fundamental mode. omodel − Flag indicating whether or not to output the model harmonic series evaluated at the times of observation in the light curve. modeloutdir − This parameter is required if the omodel flag is set. The model will be output to the file $modeloutdir/$basename.injectharm.model where $basename is the base filename of the light curve (with directory names stripped).
1 − I n j e c t t r a n s i t <" P l i s t " [ " column " c o l ] | " P f i x " p e r | " Pexpr " e x p r | " Prand " minp maxp 3 | " P l o g r a n d " minp maxp | " r a n d f r e q " minf maxf | " l o g r a n d f r e q " minf maxf> 5 <" R p l i s t " [ " column " c o l ] | " R p f i x " Rp | " Rpexpr " e x p r | " Rprand " minRp maxRp | " Rplogrand " minRp maxRp>  − Take the period of the injected transit from the input light curve list file, optionally specifying the column number. "Pfix" − Fix the period of the injected transit to the specified value. "Pexpr" − Evaluate an analytic expression for each light curve to determine the period. "Prand" − Set the period of the inject transit to a random number drawn from a uniform distribution between minp and maxp. "Plogrand" − Choose a random number for the period from a distribution that is uniform in logarithm between minp and maxp. "randfreq" − Choose a random number for the period from a distribution that is uniform in frequency between minf and maxf (in cycles per day). "lograndfreq" − Choose a random number for the period from a distribution that is uniform in the logarithm of the frequency between minf and maxf (in cycles per day). "Rplist" − Take the radius of the planet (in Jupiter radius units) from the input light curve list file, optionally specifying the column number. "Rpfix" − Fix the radius of the planet to the specified value. "Rpexpr" − Evaluate an analytic expression for each light curve to determine the planet radius. "Rprand" − Set the radius of the planet to a random number drawn from a uniform distribution between minRp and maxRp. "Rplogrand" − Set the radius of the planet to a random number drawn from a uniform-log distribution between minRp and maxRp. "Mplist" − Take the mass of the planet (in Jupiter mass units) from the input light curve list file, optionally specifying the column number. "Mpfix" − Fix the mass of the planet to the specified value. "Pexpr" − Evaluate an analytic expression for each light curve to determine the planet mass. "Mprand" − Set the mass of the planet to a random number drawn from a uniform distribution between minMp and maxMp. "Mplogrand" − Set the mass of the planet to a random number drawn from a uniform-log distribution between minMp and maxMp. "phaselist" − Take the phase of the orbit (phase= 0 corresponds to transit center) at time T = 0 from the input light curve list file, optionally specifying the column number. "phasefix" − Fix the phase at T = 0 to the specified value. "phaseexpr" − Evaluate an analytic expression for each light curve to determine the phase. "phaserand" − Set the phase at T = 0 to a random number between 0 and 1. "sinilist" − Take the sin(i), where i is the inclination of the orbit, from the input light curve list file, optionally specifying the column number. "sinifix" − Fix sin(i) to the specified value. "siniexpr" − Evaluate an analytic expression for each light curve to determine sin(i). "sinirand" − Choose a random number for sin(i) drawn from a uniform orientation distribution with the constraint that there must be a transit. "eomega" − Keyword to indicate that the eccentricity and argument of periastron will be specified. "elist" − Take the eccentricity from the input light curve list file, optionally specifying the column number. "efix" − Fix the eccentricity to the specified value. "eexpr" − Evaluate an analytic expression for each light curve to determine the eccentricity. "erand" − Adopt a random number for the eccentricity. "olist" − Take the argument of periastron (in degrees) from the input light curve list file, optionally specifying the column number. "ofix" − Fix the argument of periastron to the specified value. "oexpr" − Evaluate an analytic expression for each light curve to determine the argument of periastron. "orand" − Adopt a random number for the argument of periastron. "hk" − Keyword to indicate that e sin ω ≡ h and e cos ω ≡ k will be specified rather than the eccentricity and argument of periastron. "hlist" − Take e sin ω from the input light curve list file, optionally specifying the column number. "hfix" − Fix e sin ω to the specified value. "hexpr" − Evaluate an analytic expression for each light curve to determine e sin ω. "hrand" − Adopt a random number for e sin ω. "klist" − Take e cos ω from the input light curve list file, optionally specifying the column number. "kfix" − Fix e cos ω to the specified value. "kexpr" − Evaluate an analytic expression for each light curve to determine e cos ω. "krand" − Adopt a random number for e cos ω. "Mstarlist" − Take the mass of the star (in solar masses) from the input light curve list file, optionally specifying the column number. "Mstarfix" − Fix the mass of the star to the specified value. "Mstarexpr" − Evaluate an analytic expression for each light curve to determine the stellar mass. "Rstarlist" − Take the radius of the star (in solar radii) from the input light curve list file, optionally specifying the column number. "Rstarfix" − Fix the radius of the star to the specified value. "Rstarexpr" − Evaluate an analytic expression for each light curve to determine the stellar period. "quad" | "nonlin" − Keywords to indicate whether a 2-parameter quadratic limb darkening law or a 4-parameter nonlinear limb darkening law should be used. "ldlist" − Take the limb darkening coefficients from the input light curve list file, optionally specifying the column number of the first coefficient (other coefficients must be in the subsequent columns). "ldfix" − Fix the limb darkening coefficients to the specified values. Two values must be supplied if the "quad" keyword was used, or four values if the "nonlin" keyword was used. "ldexpr" − Evaluate an analytic expression to determine the limb darkening coefficients. Two expressions must be supplied if the "quad" keyword was used, or four expressions if the "nonlin" keyword was used. "dilute" − An optional keyword indicating that the transit signal should be scaled by a value between 0 and 1. One must then either use "list" to indicate that the scaling factor should be read from the input light curve list, "fix" to specify on the command line the value to use for all light curves, or "expr" to provide an analytic expression to be evaluated for each light curve. "omodel" − An optional keyword used to output the model transit light curve, evaluated at the times of observation in the input light curve. The model will be written to a file of the name $modeloutdir/$basename.injecttransit.model where $basename is the base filename of the input light curve, stripped of any directories.
− J s te t t i m e s c a l e d a t e s Table A.25: Input Parameters for "-Jstet" Command (Section 2.2.2.5) timescale − The time, in minutes (assuming input times are in days), that distinguishes between "near" and "far" observations. dates − A file containing a list of all times that appear in the light curves in the first column. This is used to calculate the maximum possible weight for a light curve. 1 −Killharm <"aov " | " l s " | " both " | " i n j e c t h a r m " | " f i x " Nper p e r 1 . . . perN

5
[ " outampphase " | " outampradphase" | " outRphi " | " outRradphi " ] [ " c l i p " v a l ] Table A.26: Input Parameters for "-Killharm" Command (Section 2.2.3.2) "aov" − Take the period for the harmonic series from the most recently executed "-aov" or "-aov harm" command. "ls" − Take the period for the harmonic series from the most recently executed "-LS" command. "both" − Use two periods, taking one from the most recently executed "-aov" or "-aov harm" command and the other from the most recently executed "-LS" command. "injectharm" − Take the period for the harmonic series from the most recently executed "-Injectharm" command. "fix" − Fix the period(s) to the specified value(s). Nper is the number of periods to use, and per1 through perN is a space-delimited list of the Nper periods. "list" − Take the period(s) from the input light curve list. Nper is the number of periods to use, and they must be in consecutive columns within the list file. You may optionally specify the column for the first period using the "column" keyword. Nharm − The number of higher-harmonics to use in the harmonic series that is fit for each period (frequencies of 2f 0 , 3f 0 , up to (N harm + 1)f 0 will be used for these harmonics, where f 0 is the fundamental frequency). Set this to zero to fit only a sinusoid. Nsubharm − The number of sub-harmonics to use in the harmonic series that is fit for each period (frequencies of f 0 /2, f 0 /3 up to f 0 /(N subharm + 1)). Set this to zero to fit only a sinusoid. omodel − A flag indicating whether or not to output the model harmonic series light curves. model outdir − This parameter is required if omodel is set to 1. The model harmonic series light curves will be written to files named $model outdir/$basename.killharm.model where $basename is the base filename of the input light curve (stripped of directories). "fitonly" − This optional keyword, if used, will cause Vartools to not subtract the harmonic series from the light curve, and instead only perform the fit reporting the coefficients. By default the harmonic series will be subtracted from the light curve before passing it on to the next command. "outampphase" − By default the a k and b k coefficients of the sin and cos functions in the harmonic series for harmonic number k are output. If the "outampphase" keyword is given, then the amplitudes of the harmonics A k = a 2 k + b 2 k and phases between 0 and 1 (φ k = atan2(−b k , a k )/2π) will be output.
"outampradphase" − If this keyword is given, then the amplitudes of the harmonics A k = a 2 k + b 2 k and phases in radians φ k = atan2(−b k , a k ) will be output. "outRphi" − If this keyword is given, then the amplitudes of the harmonics relative to the fundamental R k1 = A k /A 1 and relative phases φ k1 = φ k − kφ 1 (between 0 and 1) will be output. Note for subharmonics k = 1/2, 1/3, etc. For the fundamental mode the amplitude A 1 and phase φ 1 will be given. "outRradphi" − Similar to "outRphi", in this case phases are in radians. "clip" − If this optional keyword is given then the model will be fit, σ-clipping with a clipping factor given by val will be applied to the residuals, and the model will be refit to the the points which passed the clipping.  function − Analytic expression to fit to the light curve magnitudes. This expression should be linear in the free parameters. An example is 'a*t∧2+b*t+c' to fit a quadratic function in time to the light curve, with 'a', 'b' and 'c' being free parameters. paramlist − A comma-delimited list of variables which are the free parameters to optimize in the function. These must all enter linearly into the function. For the previous example one would supply 'a,b,c' for this parameter. "modelvar" − Optional keyword used to store the best-fit model light curve in a variable indicated by varname. "correctlc" − If this keyword is given then the best-fit model will be subtracted from the light curve before passing it on to the next command. By default the model is not subtracted from the light curve. "omodel" − Optional keyword used to indicate that the best-fit model light curve should be output to a file. The name of the file is given by $model outdir/$basename.linfit.model where $basename is the base filename of the input light curve stripped of directories. Theformatkeyword may be used to change the rule for naming the file as indicated below. "format" − Optional keyword used in conjunction with the "omodel" keyword to change the rule for naming the output model light curve files. In this case the file will be named $model outdir/$nameformat where instances of %s in $nameformat are replaced by $basename, instances of %d are replaced by the light curve number (starting with 1), instances of %0nd where n is an integer are replaced with the formatted light curve number, and instances of %% are replaced with %.  The maximum period to search (in the same time units as the input light curve). subsample − The periodogram will be scanned at a frequency resolution of ∆f = subsample/T where T is the time base-line of the input light curve. operiodogram − Flag used to indicate whether or not the periodograms should be output to separate files. outdir − This parameter is supplied only if operiodogram is set to 1. The periodograms will be output to files named $outdir/$basename.ls where $basename is the base filename of the input light curve stripped of any directories. "noGLS" − If this keyword is given then the traditional L-S periodogram will be calculated. By default the Generalized L-S periodogram is calculated. "whiten" − Keyword used to whiten the light curve and recalculate the periodogram for each peak. By default this is not done. "clip" − Adjust the clipping performed on the periodogram in calculating the spectroscopic S/N of a peak. By default iterative 5σ clipping is performed. clip − The sigma-clipping factor to use. clipiter − Flag which if set causes iterative clipping to be performed. "fixperiodSNR" − Report the periodogram value, the false alarm probability and the S/N for a fixed period in addition to the highest peaks found. "aov" − Use the peak period from the last executed -aov or -aov harm command for the fixperiodSNR. "ls" − Use the peak period from the last executed -LS command. "injectharm" − Use the period from the last executed -Injectharm command (the first one is used if multiple periods were injected).  Table A.29: Input Parameters for "-MandelAgolTransit" Command (Section 2.2.3.4) "bls" − Initialize the parameters based on the highest peak found with the most recent prior "-BLS" command. "blsfixper" − Initialize the parameters based on the highest peak found with the most recent prior "-BLSFixPer" command. P0 − Initial orbital period. T00 − Initial transit epoch (time at the center of a transit). r0 − Initial ratio of the planet radius to the star radius (R P /R⋆). a0 − Initial ratio of the semi-major axis to the stellar radius (a/R⋆). "i" − Keyword used to indicate that the next parameter given is the initial orbital inclination angle, in degrees. "b" − Keyword used to indicate that the next parameter given is the initial normalized impact parameter. This is the distance between the projected centers of the planet and star at conjunction divided by the sum of the planet and stellar radii (note, not just the stellar radius). "e0" − Initial orbital eccentricity. "omega0" − Initial argument of periastron in degrees. "mconst0" − Initial out-of-transit magnitude. If a negative value is specified, then the optimal value will be determined automatically. "quad" − Keyword used to indicate that a two-parameter quadratic limb darkening law is used. "nonlin" − Keyword used to indicate that a four-parameter non-linear limb darkening law is used. The model will be written to the file given by $model outdir/$basename.mandelagoltransit.model where $basename is the base filename of the input light curve, stripped of directories. This parameter should be given if and only if omodel is set to 1. "modelvar" − Optional keyword used to store the best-fit model light curve in a variable indicated by var. "ophcurve" − A keyword used to output best-fit model light curves that are uniformly sampled in phase. curve outdir − The models will be written out to files named $outdir/$basename.mandelagoltransit.phcurve. phmin − The starting phase to use in the model. phmax − The ending phase to use in the model. phstep − The phase step-size to use in the model. "ojdcurve" − A keyword used to output best-fit model light curves that are uniformly sampled in time between the first and last observed times in the light curve. curve outdir − The models will be written out to files named $outdir/$basename.mandelagoltransit.jdcurve. jdstep − The time step-size to use in the model.
− m e d i a n f i l t e r time [ " a v e r a g e " | " w e i g h t e d a v e r a g e " ] [ " r e p l a c e " ] The time separation used by the filter. By default, a high-pass filter is applied, whereby the median magnitude of all points within 'time' of a given observation is subtracted from the magnitude for that observation. "average" − Optional keyword used to take the average of the magnitudes, rather than the median. "weightedaverage" − Optional keyword used to take the weighted average of the magnitudes, rather than the median. "replace" − Optional keyword used to perform a low-pass filter, whereby the magnitude of a given observation will be replaced by the median (or average, or weighted average) magnitude of all points within 'time' of that observation.  Table A.31: Input Parameters for "-microlens" Command (Section 2.2.3.5) "f0" − Keyword used to change the method for initializing and fitting the f 0 parameter in equation 64. If this keyword is not used, the parameter will be automatically initialized and then varied in the fit. "fix" − Keyword indicating that the initial value should be fixed to the value supplied for fixval. "list" − Keyword indicating that the initial value should be taken from the input light curve list file (the column to read it from may also be optionally specified). "fixcolumn" − Keyword indicating that the initial value should be taken from the output of a previously executed command. "auto" − Keyword indicating that the initial value should be determined automatically (the default behavior). "step" − Keyword used to change the initial step-size for varyig the parameter. By default the step-size is determined automatically. "novary" − If this keyword is specified than the parameter will not be varied in the fit. "f1" − Keyword used to change the method for initializing and fitting the f 1 parameter in equation 64. If this keyword is not used, the parameter will be automatically initialized and then varied in the fit. The syntax following this keyword is the same as for "f0". "u0" − Keyword used to change the method for initializing and fitting the u 0 parameter in equation 66. If this keyword is not used, the parameter will be automatically initialized and then varied in the fit. The syntax following this keyword is the same as for "f0". "t0" − Keyword used to change the method for initializing and fitting the t 0 parameter in equation 66. If this keyword is not used, the parameter will be automatically initialized and then varied in the fit. The syntax following this keyword is the same as for "f0". "tmax" − Keyword used to change the method for initializing and fitting the tmax parameter in equation 66. If this keyword is not used, the parameter will be automatically initialized and then varied in the fit. The syntax following this keyword is the same as for "f0". "correctlc" − If this keyword is given then the best-fit model will be subtracted from the light curve before passing it on to the next command. By default the light curve is not modified by this command. "omodel" − If this keyword is given the the best-fit model, evaluated at the observed times in the light curve, will be output to a file whose name is given by $outdir/$basename.microlens where $basename is the base filename of the light curve, stripped of directories.  Table A.32: Input Parameters for "-nonlinfit" Command (Section 2.2.3.6) function − Analytic expression to fit to the light curve magnitudes. For example this might be 'a*exp(-(t-t0)∧2/2/sigma∧2)+b' for a Gaussian function in time with free parameters 'a', 't0', 'sigma' and 'b'. paramlist − A comma-delimited list of parameters to vary in the fit. For each parameter you must specify the initial guess and step-size using the format "var=init:step". For example, if the initial value for 't0' is 5.0, and its step-size is 2.0, and if the initial value for 'sigma' is 10.0 and its step-size is 8.0, you would use the expression 't0=5.0:2.0,sigma=10.0:8.0'. The initial values and step-sizes may also be more complicated analytic expressions using variables, functions, etc. Note that the free parameters should have names that are not used by any vector variables (e.g., 't', 'mag', 'err', or other variables defined by the -expr command or -inputlcformat option). Note also that these variables may be used by other commands as well (e.g., on the right-hand-side of the expression provided to the -expr command). "linfit" − Optional keyword used to indicate nay parameters in the function which enter linearly and should be optimized using linear least squares. In this case linfitparams is a comma-separated list of the linear parameters. For example, in the Gaussian case this could be 'linear a,b'. "errors" − Optional keyword used to change the magnitude uncertainties used in calculating the likelihood function (Eq. 67). error expr is an analytic expression. Parameters which are varied in the fit may also be included in this expression. "covariance" − Optional keyword used to indicated correlated uncertainties (by default errors are assume to be uncorrelated). "squareexp" − Use a square-exponential covariance matrix (Eq. 69). amp var − The name of the Vartools variable which should be used as a in equation 69. If this variable does not appear in paramlist then it should be specified here as 'amp var=expr' where amp var is the name of the variable to use, and expr is an expression used to determine the fixed value to be used for this parameter. rho var − The name of the Vartools variable which should be used as ρ in equation 69. If this variable does not appear in paramlist then it should be specified here as 'rho var=expr' where rho var is the name of the variable to use, and expr is an expression used to determine the fixed value to be used for this parameter. "exp" − Use a exponential covariance matrix (Eq. 70). amp var − The name of the Vartools variable which should be used as a in equation 70. If this variable does not appear in paramlist then it should be specified here as 'amp var=expr' where amp var is the name of the variable to use, and expr is an expression used to determine the fixed value to be used for this parameter. rho var − The name of the Vartools variable which should be used as ρ in equation 70. If this variable does not appear in paramlist then it should be specified here as 'rho var=expr' where rho var is the name of the variable to use, and expr is an expression used to determine the fixed value to be used for this parameter. "matern" − Use a Matérn covariance matrix (Eq. 71). amp var − The name of the Vartools variable which should be used as a in equation 71. If this variable does not appear in paramlist then it should be specified here as 'amp var=expr' where amp var is the name of the variable to use, and expr is an expression used to determine the fixed value to be used for this parameter.
rho var − The name of the Vartools variable which should be used as ρ in equation 71. If this variable does not appear in paramlist then it should be specified here as 'rho var=expr' where rho var is the name of the variable to use, and expr is an expression used to determine the fixed value to be used for this parameter. nu var − The name of the Vartools variable which should be used as ν in equation 71. If this variable does not appear in paramlist then it should be specified here as 'nu var=expr' where nu var is the name of the variable to use, and expr is an expression used to determine the fixed value to be used for this parameter. "priors" − Keyword used to indicate a list of priors to place on the variables. priorlist is a comma-separated list, where each entry in the list is an analytic expresion that should evaluate to −2 ln(P ) where P is the prior probability for a variable given its value. For example, to place a Gaussian prior on 't0' with mean 4.0 and standard deviation 3.0 you would use 'prior (t0-4.0)∧2/3.0∧2'. "constraints" − Use this keyword to place constraints on the parameters. constraintlist is a comma-separated list of expressions. For example, to require a positive 'sigma' you can use 'constraints sigma>0'. "amoeba" − Use the down-hill simplex optimization algorithm. "tolerance" − Optional keyword to change the convergence tolerance. Here 'tol' corresponds to TOL in equation 72. "maxsteps" − Optional keyword to change the maximum number of amoeba iterations that are tried before giving up. "mcmc" − Run a Differential Evolution Markov Chain Monte Carlo procedure. "Naccept" − Optionally specify the number of accepted links to run in a given fit. Once this number is reached, the MCMC will terminate. "Nlinkstotal" − Optionally specify a maximum total number of links to run. Once this number is reached, the MCMC will terminal. This option is assumed by default, with a limit of 100,000 links. "fracburnin" − Optionally change the initial fraction of the chain that is ignored in computing statistics to be reported from the posterior distribution (the default is 0.1). "eps" − Optionally change the value of ǫ to use in equation 73. The default value is 0.001. "skipamoeba" − By default the down-hill simplex algirhtm is run initially before starting the mcmc procedure. This can be skipped by using this keyword. "chainstats" − By default the median and standard deviation for each of the parameters varied in the MCMC will be included in the output table. You may use this keyword to change the statistics, and/or the quantities that are used. Here 'exprlist' is a comma-separated list of analytic expressions to calculate from the chain, and 'statslist' is a comma-separated list of statistics to report for each of the expressions. The available statistics are the same as for the -stats command. "maxmemstore" − Change the maximum limit on the total amount of memory to be used by this MCMC command. This will limit the length of the chain used for calculating the output statistics. 'maxmem' is the limit in GB. The default is 4.0. "outchains" − Use this keyword to output the MCMC chains. The chains will be output to files named $model outdir/$basename.mcmc, where $basename is the base filename of the input light curve file, stripped of directories. "format" − Use this keyword to change the naming convention for the output MCMC files. In this case the file will be named $model outdir/$nameformat where instances of %s in $nameformat are replaced by $basename, instances of %d are replaced by the light curve number (starting with 1), instances of %0nd, where n is an integer, are replaced with the formatted light curve number, and instances of %% are replaced with %. "printevery" − Use this keyword to change the number of links that are written to the output file. By default every link is written out. If this keyword is used, then every 'N'th link is written out. "modelvar" − Optional keyword used to store the best-fit model light curve in a variable indicated by varname. "correctlc" − Use this keyword to subtract the best-fit model light curve from the input light curve before passing it on to the next command. "omodel" − If this keyword is given, the best-fit model, evaluated at the observed times in the light curve, will be output to a file whose name is given by $model outdir/$basename.nonlinfit.model where $basename is the base filename of the light curve, stripped of directories. "format" − Use this keyword to change the naming convention for the output model light curve. In this case the file will be named $model outdir/$nameformat where instances of %s in $nameformat are replaced by $basename, instances of %d are replaced by the light curve number (starting with 1), instances of %0nd, where n is an integer, are replaced with the formatted light curve number, and instances of %% are replaced with %.
outdir − The directory to which the light curve files will be written. By default the output filenames will be $outdir/$basename where $outdir is the value supplied for 'outdir' and $basename is the base filename of the light curve, stripped of directories. The supplied parameter will be interpreted as the output directory only if a light curve list was used for input (-l). outname − The name of the output light curve file. This will be used if an individual light curve file was used for input (-i). "nameformat" − Use this keyword to change the naming convention for the output light curves. In this case the file will be named $outdir/$formatstring where instances of %s in $formatstring are replaced by $basename, instances of %d are replaced by the light curve number (starting with 1), instances of %0nd, where n is an integer, are replaced with the formatted light curve number, and instances of %% are replaced with %. "columnformat" − Change the format of the output light curves. By default they will have three columns: time, mag and err. Here 'formatstring' is a comma-separated list of variable names to output, optionally with a ':' after each variable name to specify the "printf"-style format to use for that variable. For example, 'columnformat t:%.17g,mag:%.5f,err:%.5f,xpos:%.3f' would output the variables t, mag, err, and xpos using 17-digit-precision for t, and floating points for mag, err and xpos, with 5, 5, and 3 digits after the decimal, respectively. Here 'xpos' is a non-default variable that one might have read-in with the -inputlcformat option, or produced through one of the other commands.
If the light curves are output in fits format, then terms after the ':' will be used to specify the units of the column to list in the light curve header. "fits" − Use this keyword to output the light curves in binary FITS table format. The output light curves will have the extension ".fits" appended to their filenames, if not already present. "noclobber" − Use this keyword to prevent over-writing any existing files. Vartools will terminate if it encounters an existing file when this keyword is used.
[ " p h a s e v a r " v a r ] [ " s t a r t p h a s e " s t a r t p h a s e ] Table A.34: Input Parameters for "-Phase" Command (Section 2.2.6.8) "aov" − Use this keyword to take the period at which to phase-fold the light curve from the most recently issued -aov or -aov harm command. "ls" − Use this keyword to take the period at which to phase-fold the light curve from the most recently issued -LS command. "bls" − Use this keyword to take the period at which to phase-fold the light curve from the most recently issued -BLS command. "T0" − By default the light curves are folded such that phase zero corresponds to time t = 0. Use this keyword to specify a different reference time for phase zero. "bls" − Take the reference time to be the time of central transit found in the most recent -BLS command.
Here 'phaseTc' is the phase to adopt for the central transit time. "phasevar" − Use this keyword to write the phases to the variable var rather than overwriting the time variable. If the variable var does not exist it will be created. If it does exist, it must be a double-precision light curve vector (e.g., it cannot be the name of an output column, or a parameter used in a fit, etc.). "startphase" − Use this keyword to change the range over which the phases are expressed. By default they run from 0 to 1. With this keyword they will run from startphase to startphase + 1.  − Perform linear interpolation between points. "spline" − Perform cubic spline interpolation. "left" − Specify the left boundary condition for the spline. "right" − Specify the right boundary condition for the spline. "splinemonotonic" − Perform cubic spline interpolation with the interpolating function constrained to be monotonic between input observations. "bspline" − Interpolate with a Basis-spline function. "nbreaks" − The number of breaks to use in the B-spline (the default is 15). If the number given is < 2, then the routine will increase the number of breaks until χ 2 per degree of freedom is less than one. Caution, this can be quite slow. "order" − The order of the spline function to use (the default is 3). "file" − Optional keyword to read the time-base on which to resample the light curve from a file. By default Vartools will resample the light curve onto a uniform time-base with the starting and stopping times equal to the minimum and maximum observed times in the light curve, respectively, and with a time step equal to the minimum time separation between consecutive points in the input light curve. "fix" − Use the same time-base file for all light curves. This file is given by 'times file'. Use the optional "column" keyword to specify the column in the file which stores the times to use. "list" − Allow a different time-base file for each light curve. The filenames in this case are taken from the input light curve list file. Use the optional "listcolumn" keyword to change the column in the input light curve list file to use, use the "tcolumn" keyword to change the columns in the files to use for the new times. "tstart" − Optional keyword to change the starting time for a regularly spaced grid of times to resample the light curve at. "expr" − In addition to the standard "fix", "fixcolumn", and "list" options for determining the starting time, one may also use this keyword to set it equal to the result from evaluating an analytic expression given by 'expression'. "tstop" − Optional keyword to change the starting time for a regularly spaced grid of times to resample the light curve at. Options for determining this parameter are the same as for "tstart". "delt" − Optional keyword to change the time step for a regularly spaced grid of times to resample the light curve at. Options for determining this parameter are the same as for "tstart". "Npoints" − Optional keyword to change the number of points to resample the light curve onto. Either "delt" or "Npoints" may be specified. Options for determining this parameter are the same as for "tstart". "gaps" − By default the same resampling method will be used for all points in the light curve. By giving the "gaps" keyword you can make the method depend on how far away a resampled time is from the closest observed time. You first need to indicate how the time separation used to distinguish between the near and far points will be determined, you then indicate the method for resampling the "far" points (the near points will be resampled using the method already indicated before providing the "gaps" keyword). Below are the options for determining the time separation. The options for determining the resampling method are the same as discussed above.
"fix" − Fix the time separation to value specified by 'time sep'. "fixcolumn" − Set the time separation equal to a value calculated by a previously executed command. "list" − Take the value from the input light curve list file. "expr" − Set the value equal to the result from evaluating an analytic expression. "frac min sep" − Set the time separation to a fixed factor times the minimum separation between subsequent points in the input light curve (e.g., if you give "frac min sep 5.0" and the minimum separation between points in the light curve is 1 day, then the separation times scale will be set to 5 days). "frac med sep" − Same as "frac min sep", except using the median separation between subsequent points rather than the minimum separation. "percentile sep" − Set the time separation to a specified percentile of the separations in the input light curve (e.g., if you give "percentile sep 70", then the N separations in the input light curve will be ordered, and the separation to be used for distinguishing between near and far points will be set equal to the 0.7 × N longest separation in the input light curve). "extrap" − Change the method for evaluating resampled points which are extrapolations rather than interpolations. Options are as already discussed. By default extrapolations and interpolations are performed using the same method.    − If this optional keyword is given, then the times specified below will be removed from the light curve. By default times specified are included, and other times are excluded. "JDrange" − Provide a fixed range of times to keep (or exclude) with 'minJD' and 'maxJD'. "JDrangebylc" − Provide a range of times to keep (or exclude), allowing the limits on the range to change between light curves. If this keyword is given, then the method for determining the 'minJD' should be provided (either fixed to a constant value, read from the input light curve list, set to the result of a previously executed command, or set equal to the result from evaluating an analytic expression. The method for determining 'maxJD' should then be provided, with the same set of options.    − Use this keyword to take the period for the harmonic series from the most recent -aov oraov harm command. "ls" − Use this keyword to take the period for the harmonic series from the most recent -LS command. "bls" − Use this keyword to take the period for the harmonic series from the most recent -BLS command. "list" − Use this keyword to take the period for the harmonic series from the input light curve list file (optionally specify the column to use).  − Take the rotation period for the model from the most recent -aov or -aov harm command. "ls" − Take the rotation period for the model from the most recent -LS command. "list" − Take the rotation period for the model from the input light curve list (optionally specifying the column to use). "fix" − Fix the rotation period for the model to the value specified by 'period'. "fixcolumn" − Set the rotation period for the model to the value from a previously executed command. a0 − The initial value to use for a in equation 76. b0 − The initial value to use for b in equation 76. alpha0 − The initial value to use for the spot angular radius in degrees. i0 − The initial value to use for the inclination of the stellar rotation axis in degrees (90 • corresponds to the rotation axis being perpendicular to the line-of-sight). chi0 − The initial value to use for the spot latitude in degrees (0 • for a spot at the equator). psi00 − The initial value to use for the longitude of the spot center (at the first time instance in the light curve) in degrees.   The number of trend vectors to remove for which the color terms are to be specified initially. These will be read-in from the input light curve list file, by default starting from the next unused column. The terms must be in Ninput color consecutive columns. You can change the column to use for the first term with the "column" keyword. Ninput airmass − The number of trend vectors to remove for which the airmass terms are to be specified initially. initial airmass file − A file with the initial airmass trends to use. The first column in this file should be the JDs (or image IDs if the -matchstringid option was given), and the subsequent Ninput airmass columns are the initial airmass trends. sigma clip1 − σ-clipping factor used in calculating the mean magnitudes of the light curves. sigma clip2 − σ-clipping factor used in determining whether or not points contribute to the airmass or color terms when performing the fit. saturation − Any points with magnitude less than 'saturation' will not contribute to the fit. correctlc − Flag indicating whether or not the best-fit trend model should be subtracted from the light curve before passing it on to the next command. omodel − Flag indicating whether or not the model light curve files, evaluated at the observed times in the input light curve, should be output.   − If this optional keyword is used, the the light curve will be simultaneously decorrelated against additional light-curve-specific signals. iterativeflag − Flag indicating whether the decorrelation and TFA will be done iteratively (flag set to 1; this is faster) or if they will be done simultaneously (flag set to 0; slower but more correct). Nlcterms − The number of decorrelation trends to use. lccolumn1 − The column in the input light curve to use for the first light-curve-specific trend to decorrelate against. lcorder1 − The polynomial order to use for the first trend. . . .
− The columns and orders should then be provided for the rest of the Nlcterms trends. pixelsep − Trend stars within pixelsep of the light curve in question will not be used in detrending the light curve. "xycol" − Optional keyword to change the columns in the input light curve list file used to read-in the X and Y coordinates of the light curves. By default these are taken from the next unused columns. correctlc − Flag indicating whether or not the best-fit trend model should be subtracted from the light curve before passing it on to the next command. ocoeff − Flag indicating whether or not the trend coefficients should be output for each light curve. coeff outdir − The trend coefficients will be written to files named $coeff outdir/$basename.tfa.coeff where $basename is the base filename of the input light curve, stripped of directories. This parameter should be given if and only if ocoeff is set to 1. omodel − Flag indicating whether or not to output the best-fit model trends, evaluated at the times in the input light curves. model outdir − The trend models will be written to files named $model outdir/$basename.tfa.model where $basename is the base filename of the input light curve, stripped of directories. This parameter should be given if and only if omodel is set to 1. dotfafirst − Flag indicating whether TFA should be applied to the input light curve first, with the signal to preserve determined on the residual in each iteration (set the flag to 1), or if the signal is determined and subtracted from the light curve first and then TFA is applied to the residual in each iteration (set the flag to 0). tfathresh − The iterations will stop if the fractional change in the r.m.s. is less than 'tfathresh'. maxiter − The iterations will stop once 'maxiter' iterations have completed. "bin" − Use the binned light curve for the model signal to preserve. nbins − The number of bins to use. "period" − Optional keyword to do phase-binning, rather than binning in time. "aov" − Take the period for phase-binning from the most recently executed -aov of -aov harm command. "ls" − Take the period for phase-binning from the most recently executed -LS command. "bls" − Take the period for phase-binning from the most recently executed -BLS command. "list" − Take the period for phase-binning from the input light curve list file (optionally specifying the column). "fix" − Fix the period for phase-binning to the value specified by 'period'. "signal" − Use a fixed signal form read-in from a file. Use 'filename' to specify a file listing the signal files, one for each light curve to process. Each signal file should contain the signal in the second column. The quantity a × S + b is fit to the light curve simultaneously with TFA, here a and b are free parameters and S is the signal. "harm" − Use a harmonic series to model the signal. This will be fit simultaneously to the light curve with TFA. Nharm − The number of harmonics to use. Nsubharm − The number of sub-harmonics to use. "period" − Use this keyword to change the period for the harmonic series (the default is the time spanned by the light curve). "aov" − Use this keyword to take the period from the most recent -aov or -aov harm command. "ls" − Use this keyword to take the period from the most recent -LS command. "bls" − Use this keyword to take the period from the most recent -BLS command. "list" − Read the period from the light curve list file (optionally specify the column). "fix" − Fix the period to a specified value.
1 −wwz <"maxfreq " <"auto " | maxfreq>> <"f r e q s a m p " freqsamp> <"tau0 " <"auto " | tau0 >> <"tau1 " <"auto " | tau1 >>  − Indicate the maximum frequency for which the wavelet transform will be calculated. "auto" − If this keyword is given then the maximum frequency will be 1/(2∆t min ) where ∆t min is the minimum time seperation between consecutive points in the light curve. maxfreq − The maximum frequency to use. "freqsamp" − Specify the frequency sampling as freqsamp/T where T is the time baseline of the light curve. "tau0" − Indicate the minimum time-shift to test (the time-shift is τ in Eq. 17). Use the "auto" keyword to set this to the minimum time in the light curve, otherwise provide a value. "tau1" − Indicate the maximum time-shift to test. Use the "auto" keyword to set this to the maximum time in the light curve, otherwise provide a value. "dtau" − Indicate the time-shift step to use. Use the "auto" keyword to set this to ∆t min , otherwise provide a value. "c" − Optionally provide a value for the c parameter in equation 18. The default value is c = (8π 2 ) −1 . "outfulltransform" − Use this keyword to output the full wavelet transform for each light curve (i.e., calculated at every trial time-shift and frequency). outdir − The directory to output the transforms to. The default naming convention is $outdir/$basename.wwz where $basename is the base filename of the input light curve, stripped of directories. "fits" − Optionally output the transforms as multi-extension FITS image files. "pm3d" − Optionally output the transforms as ascii tables in a format suitable for plotting with the gnuplot pm3d plotting style. "format" − Modify the naming convention for the files. In this case the file will be named $outdir/$format where instances of %s in $format are replaced by $basename, instances of %d are replaced by the light curve number (starting with 1), instances of %0nd, where n is an integer, are replaced with the formatted light curve number, and instances of %% are replaced with %. "outmaxtransform" − Use this keyword to output the transform that maximizes Z (Eq. 29) over frequencies as a function of time-shift. outdir − The directory to output these files to. The default naming convention is $outdir/$basename.mwwz where $basename is the base filename of the input light curve, stripped of directories. "format" − Modify the naming convention for the files. In this case the file will be named $outdir/$format where instances of %s in $format are replaced by $basename, instances of %d are replaced by the light curve number (starting with 1), instances of %0nd, where n is an integer, are replaced with the formatted light curve number, and instances of %% are replaced with %.

Appendix B. Columns Output By Each Command
The tables below list the columns added by each command to the output table. The columns will be named $colnum" "$basecolname" "$commandnum where $colnum" " is prepended only if the -numbercolumns option is given to Vartools, $colnum is the column number (starting at 1), $basecolname is as listed below, and $commandnum is the number of the command that this column was produced by (i.e., if one runs "vartools -i input.txt -rms -chi2 -rms" the columns produced by the first "-rms" call would have $commandnum= 1, those produced by the "-chi2" command would have $commandnum= 2, and those produced by the second "-rms" call would have $commandnum= 3). When using a "fixcolumn" option to one of the Vartools commands to set a parameter equal to the value of the command, one can either provide the column number ($colnum), or the column name without the prepended column number ($basecolname" "$commandnum). Table B.46: Output Columns for "-addnoise" Command (Section 2.2.5.1) Addnoise Gamma 1 − The value used for γ, as discussed in section 2.2.5. This is only output if the "wavelet" noise model is used, and if "list" is used to initialize γ. Addnoise Sig Red 1 − The value used for the red-noise component standard deviation. This is only output if "list" is used to initialize red-noise standard deviation. Addnoise Sig White 1 − The value used for the white-noise component standard deviation.. This is only output if "list" is used to initialize the white noise standard deviation. Addnoise Rho 1 − The value used for ρ, as discussed in section 2.2.5. This is only output if the "squareexp", "exp" or "matern" noise model is used, and if "list" is used to initialize ρ. Addnoise Nu 1 − The value used for ν, as discussed in section 2.2.5. This is only output if the "matern" noise model is used, and if "list" is used to initialize ν.  The period for peak number $peaknum. $peaknum runs from 1 for the highest peak in the periodogram to 'Npeak', the number of peaks to identify as given on the command line. Note that all columns for a given peak are output together, and then the same quantities are repeated for the next peak, until all 'Npeak' peaks have been reported. AOV LOGSNR $peaknum − The spectroscopic S/N measured on the logarithmic AoV periodogram for peak $peaknum. This is output only if the "uselog" keyword was given to the -aov command. AOV $peaknum − The value of the AoV periodogram for peak $peaknum. This is output only if the "uselog" keyword was not given to the -aov command. AOV SNR $peaknum − The spectroscopic S/N measured on the AoV periodogram for peak $peaknum. This is output only if the "uselog" keyword was not given to the -aov command. AOV NEG LN FAP $peaknum − The negative natural logarithm of the formal false alarm probability for peak $peaknum. This is output only if the "uselog" keyword was not given to the -aov command. Mean lnAOV $peaknum − The mean value of the logarithmic AoV periodogram at the whitening cycle used to find peak $peaknum. This is output only if the "uselog" and the "whiten" keywords were both given to the -aov command. RMS lnAOV $peaknum − The standard deviation of the logarithmic AoV periodogram at the whitening cycle used to find peak $peaknum. This is output only if the "uselog" and the "whiten" keywords were both given to the -aov command. PeriodFix − The fixed period at which to report the periodogram properties. This is output only if the "fixperiodSNR" keyword was given to the -aov command. AOV LOGSNR PeriodFix − The spectroscopic S/N measured on the logarithmic AoV periodogram for the fixed period. This is output only if the "fixperiodSNR" and the "uselog" keywords were both given to the -aov command. AOV PeriodFix − The value of the AoV periodogram at the fixed period. This is output only if the "fixperiodSNR" keyword is given and the "uselog" keyword is not given. AOV SNR PeriodFix − The spectroscopic S/N of the AoV periodogram at the fixed period. This is output only if the "fixperiodSNR" keyword is given and the "uselog" keyword is not given.

AOV NEG LN FAP PeriodFix −
The negative natural logarithm of the formal false alarm probability at the fixed period. This is output only if the "fixperiodSNR" keyword is given and the "uselog" keyword is not given. Mean lnAOV − The mean value of the logarithmic AoV periodogram. This is output only if the "uselog" keyword is given and the "whiten" keyword is not given. RMS lnAOV − The standard deviation of the logarithmic AoV periodogram. This is output only if the "uselog" keyword is given and the "whiten" keyword is not given. Table B.49: Output Columns for "-aov harm" Command (Section 2.2.1.5) Period $peaknum − The period for peak number $peaknum. $peaknum runs from 1 for the highest peak in the periodogram to 'Npeak', the number of peaks to identify as given on the command line. Note that all columns for a given peak are output together, and then the same quantities are repeated for the next peak, until all 'Npeak' peaks have been reported. AOV HARM $peaknum − The value of the AoV periodogram for peak $peaknum. This is output only if a positive number was specified for the 'Nharm' parameter input to the -aov harm command. AOV HARM NEG LOG FAP $peaknum − The negative logarithm of the formal false alarm probability for peak $peaknum. This is output here only if a value less than or equal to zero was specified for the 'Nharm' parameter input to the -aov harm command. AOV HARM NHARM $peaknum − The optimum number of harmonics used for peak $peaknum. This is output here only if a value less than or equal to zero was specified for the 'Nharm' parameter input to the -aov harm command. AOV HARM SNR $peaknum − The spectroscopic S/N measured on the AoV periodogram for peak number $peaknum. AOV HARM NEG LOG FAP $peaknum − The negative logarithm of the formal false alarm probability for peak $peaknum. This is output here only if a positive number was specified for the 'Nharm' parameter input to the -aov harm command. Mean AOV HARM $peaknum − The mean value of the AoV periodogram for the whitening cycle in which peak $peaknum was found. This is only output if the "whiten" keyword was given. RMS AOV HARM $peaknum − The standard deviation of the AoV periodogram for the whitening cycle in which peak $peaknum was found. This is only output if the "whiten" keyword was given.

−
The fixed period at which to report the periodogram properties. This is output only if the "fixperiodSNR" keyword was given to the -aov harm command. AOV HARM PeriodFix − The value of the AoV periodogram at the fixed period. This is output only if the "fixperiodSNR" keyword is given. AOV HARM SNR PeriodFix − The spectroscopic S/N of the AoV periodogram at the fixed period. This is output only if the "fixperiodSNR" keyword is given. AOV HARM NEG LN FAP PeriodFix − The negative logarithm of the formal false alarm probability for the fixed period. This is output here only if the "fixperiodSNR" keyword is given, and if a positive value was specified for the 'Nharm' parameter input to the -aov harm command. Mean AOV HARM − The mean of the AoV periodogram. This is output only if the "whiten" keyword was not provided. RMS AOV HARM − The standard deviation of the AoV periodogram. This is output only if the "whiten" keyword was not provided. The value of ∆χ 2 = χ 2 0 − χ 2 invtr where χ 2 0 is for a model with no variation, and χ 2 invtr is the best-fit inverse transit model found.

BLS MeanMag
− The mean magnitude of the light curve which is subtracted before applying the BLS algorithm. The value of ∆χ 2 = χ 2 0 − χ 2 invtr where χ 2 0 is for a model with no variation, and χ 2 invtr is the best-fit inverse transit model found at the fixed period. BLSFixPer MeanMag − The mean magnitude of the light curve which is subtracted before calculating the BLS statistics. The fixed transit depth adopted for this command. This is output here only if the "fixdepth" keyword was used. BLSFixDurTc Qingress − The fixed fractional ingress duration (duration of ingress divided by total transit duration) adopted for this command. This is output only if the "fixdepth" keyword was used. If the "qgress" keyword was not given, this will have a value of 0. BLSFixDurTc Period $peaknum − The period for peak number $peaknum. $peaknum runs from 1 for the highest S/N peak in the periodogram to 'Npeak', the number of peaks to identify as given on the command line. Note that all columns for a given peak are output together, and then the same quantities are repeated for the next peak, until all 'Npeak' peaks have been reported.
BLSFixDurTc SN $peaknum − The spectroscopic S/N estimate for peak number $peaknum equal to the SDE (Eq. 12), except that a 3σ clipping is first applied to the BLS spectrum before computing the average value and standard deviation to use in that equation. The BLS peaks are sorted by this quantity. BLSFixDurTc SR $peaknum − The Signal Residue for peak $peaknum (Eq. 10). BLSFixDurTc SDE $peaknum − The Signal Detection Efficiency for peak $peaknum (Eq. 12). BLSFixDurTc Depth $peaknum − The transit depth in magnitudes (or whatever units are used for the input light curve) for peak number $peaknum. BLSFixDurTc Qtran $peaknum − The fractional transit duration (duration divided by orbital period) for peak number $peaknum. BLSFixDurTc Qingress $peaknum − The fractional ingress duration (duration of ingress divided by total transit duration) for peak number $peaknum. This is output only if the "fittrap" keyword was used and the "fixdepth" keyword was not used. BLSFixDurTc OOTmag $peaknum − The out-of-transit mangitude level for peak number $peaknum. This is output here only if the "fittrap" keyword was used and the "fixdepth" keyword was not used. BLSFixDurTc deltaChi2 − The value of ∆χ 2 = χ 2 0 − χ 2 tr for peak number $peaknum where χ 2 0 is for a model with no variation, and χ 2 tr is for the transit model.  The mean magnitude of the light curve, using the input uncertainties to determine the weights. This is the value used in calculating the magnitude residuals for computing χ 2 .  The frequency in cycles per day (assuming input times are in days) for peak number $peaknum identified in the "dirty" DFT spectrum (i.e., the spectrum before application of the "CLEAN" deconvolution algorithm). $peaknum runs from 1 for the highest peak in the spectrum to 'Npeaks', the number of peaks to identify in the "dirty" spectrum as given on the command line. Note that all columns for a given peak are output together, and then the same quantities are repeated for the next peak, until all 'Npeaks' peaks have been reported. These columns are only reported if the "findddirtypeaks" keyword is given.
The value of the power spectrum at peak number $peaknum. This is equal to P (f ) = F 2 r (f ) + F 2 i (f ) where Fr and F i are the real and imaginary components, respectively, of F (f ) given in equation 7. If the "useampspec" keyword is given on the command line, then it is the square root of this quantity that is reported. DFTCLEAN DSPEC PEAK SNR $peaknum − The spectroscopic S/N ratio for peak number $peaknum. This equal to (P (f ) −P )/σ P , where P (f ) is the value of the power spectrum at peak number $peaknum (or amplitude spectrum if the "useampspec" keyword is given), andP and σ P are the average and standard deviation of the spectrum, respectively, after σ-clipping. DFTCLEAN DSPEC AVESPEC − The value ofP used in computing the S/N ratio(s) for the dirty power spectrum. This is output only if the "verboseout" keyword is given, and if the "finddirtypeaks" keyword is given with 'Npeaks' greater than 0. DFTCLEAN DSPEC STDSPEC − The value of σ P used in computing the S/N ratio(s) for the dirty power spectrum. This is output only if the "verboseout" keyword is given, and if the "finddirtypeaks" keyword is given with 'Npeaks' greater than 0. DFTCLEAN DSPEC AVESPEC NOCLIP − The value ofP for the dirty spectrum when no clipping is applied. This is output only if the "verboseout" keyword is given, and if the "finddirtypeaks" keyword is given with 'Npeaks' greater than 0. DFTCLEAN DSPEC STDSPEC NOCLIP − The value of σ P for the dirty spectrum when no clipping is applied. This is output only if the "verboseout" keyword is given, and if the "finddirtypeaks" keyword is given with 'Npeaks' greater than 0. DFTCLEAN CSPEC PEAK FREQ $peaknum − The frequency in cycles per day (assuming input times are in days) for peak number $peaknum identified in the "CLEAN" DFT spectrum. $peaknum runs from 1 for the highest peak in the spectrum to 'Npeaks', the number of peaks to identify in the "CLEAN" spectrum as given on the command line. Note that all columns for a given peak are output together, and then the same quantities are repeated for the next peak, until all 'Npeaks' peaks have been reported. These columns are only reported if the CLEAN algirthm is performed, and the "findcleanpeaks" keyword is given. DFTCLEAN CSPEC PEAK POW $peaknum − The value of the "CLEAN" power spectrum at peak number $peaknum. DFTCLEAN CSPEC PEAK SNR $peaknum − The spectroscopic S/N ratio for peak number $peaknum in the "CLEAN" spectrum. DFTCLEAN CSPEC AVESPEC − The value ofP used in computing the S/N ratio(s) for the clean power spectrum. This is output only if the "verboseout" keyword is given, and if the "findcleanpeaks" keyword is given with 'Npeaks' greater than 0. DFTCLEAN CSPEC STDSPEC − The value of σ P used in computing the S/N ratio(s) for the clean power spectrum. This is output only if the "verboseout" keyword is given, and if the "findcleanpeaks" keyword is given with 'Npeaks' greater than 0. DFTCLEAN CSPEC AVESPEC NOCLIP − The value ofP for the clean spectrum when no clipping is applied. This is output only if the "verboseout" keyword is given, and if the "findcleanpeaks" keyword is given with 'Npeaks' greater than 0. DFTCLEAN CSPEC STDSPEC NOCLIP − The value of σ P for the clean spectrum when no clipping is applied. This is output only if the "verboseout" keyword is given, and if the "findcleanpeaks" keyword is given with 'Npeaks' greater than 0. The value of χ 2 new /χ 2 old where χ 2 old and χ 2 new are the reduced χ 2 values of the light curve before and after transforming the magnitude uncertainties, respectively. The name of the light curve with the largest amplitude (in flux) harmonic series at the indicated period. This light curve is the most likely source for the measured variability. FindBlends FluxAmp − The amplitude (in flux) of the highest-amplitude light curve. Table B.68: Output Columns for "-fluxtomag" Command (Section 2.2.6.7) No columns added to the output table. Table B.69: Output Columns for "-GetLSAmpThresh" Command (Section 2.2.1.2) LS AmplitudeScaleFactor − The scaling factor α for which the detected FAP would be equal to the threshold. LS MinimumAmplitude − The minimum peak-to-peak amplitude that the light curve could have had, and still have been detected by L-S with a FAP below the threshold.  The orbital phase at time t = 0 (phase= 0 corresponds to transit center) for the transit signal added to the light curve. Injecttransit sin i − The sine of the orbital inclination angle for the transit signal added to the light curve. Injecttransit e − The orbital eccentricity e for the transit signal added to the light curve. This is output only if the "eomega" keyword is used. Injecttransit omega − The argument of periastron ω (in degrees) for the transit signal added to the light curve. This is output only if the "eomega" keyword is used. Injecttransit h − The value of e sin ω for the transit signal added to the light curve. This is output only if the "hk" keyword is used. Injecttransit k − The value of e cos ω for the transit signal added to the light curve. This is output only if the "hk" keyword is used. Injecttransit Mstar − The mass of the star in Solar masses for the transit signal added to the light curve. Injecttransit Rstar − The radius of the star in Solar radii for the transit signal added to the light curve. Injecttransit ld $i − The value of limb darkening coefficient $i ($i takes values from 1 to 2 for the quadratic limb-darkening law, or from 1 to 4 for the non-linear limb-darkening law) for the transit signal added to the light curve.  Here $i takes values from 1 to 'Nper', where 'Nper' is the number of periods indicated on the command-line ('Nper' is 1 for the "aov", "ls" and "injectharm" options, it is 2 for the "both" option). All columns for a given value of $i are output before the columns for the next value of $i are output. Killharm Per$i Subharm $k Sincoeff − The fitted value of c i,k from equation 49 where k =$k. $k takes values from 2 to 'Nsubharm'+1, where 'Nsubharm' is the number of subharmonics to use as specified on the command line. This is output only if none of the keywords "outampphase", "outampradphase", "outRphi" or "outRradphi" are given. Killharm Per$i Subharm $k Coscoeff − The fitted value of d i,k from equation 49 where k =$k. $k takes values from 2 to 'Nsubharm'+1, where 'Nsubharm' is the number of subharmonics to use as specified on the command line. This is output only if none of the keywords "outampphase", "outampradphase", "outRphi" or "outRradphi" are given. Killharm Per$i Subharm Amp $k − The fitted value of B i,k from equation 51. This is output only if one of the keywords "outampphase" or "outampradphase" is given. Killharm Per$i Subharm Phi $k − The fitted value of ψ i,k from equation 55. This is output only if one of the keywords "outampphase" or "outampradphase" is given. If "outampphase" is given, then ψ i,k /(2π) is reported. Killharm Per$i Subharm R $k 1 − The fitted value of the relative amplitude Q i,k,1 from equation 53. This is output only if the "outRphi" or "outRradphi" keyword is given. Killharm Per$i Subharm Phi $k 1 − The fitted value of the relative phase ψ i,k,1 from equation 57. This is output only if the "outRphi" or "outRradphi" keyword is given. If "outRphi" is given, then ψ i,k,1 /(2π) is reported. Killharm Per$i Fundamental Sincoeff − The fitted value of a i,1 from equation 49. This is output only if none of the keywords "outampphase", "outampradphase", "outRphi" or "outRradphi" are given. Killharm Per$i Fundamental Coscoeff − The fitted value of b i,1 from equation 49. This is output only if none of the keywords "outampphase", "outampradphase", "outRphi" or "outRradphi" are given. Killharm Per$i Fundamental Amp − The fitted value of A i,1 from equation 50. This is output only if one of the keywords "outampphase", "outampradphase", "outRphi" or "outRradphi" is given. Killharm Per$i Fundamental Phi − The fitted value of φ i,1 from equation 54. This is output only if one of the keywords "outampphase", "outampradphase", "outRphi" or "outRradphi" is given. If "outampphase" or "outRphi" is given, then 'φ i,1 /(2π) is reported. Killharm Per$i Harm $k Sincoeff − The fitted value of a i,k from equation 49 where k =$k. $k takes values from 2 to 'Nharm'+1, where 'Nharm' is the number of harmonics to use as specified on the command line. This is output only if none of the keywords "outampphase", "outampradphase", "outRphi" or "outRradphi" are given. Killharm Per$i Harm $k Coscoeff − The fitted value of b i,k from equation 49 where k =$k. $k takes values from 2 to 'Nharm'+1, where 'Nharm' is the number of harmonics to use as specified on the command line. This is output only if none of the keywords "outampphase", "outampradphase", "outRphi" or "outRradphi" are given. Killharm Per$i Harm Amp $k − The fitted value of A i,k from equation 50. This is output only if one of the keywords "outampphase" or "outampradphase" is given. Killharm Per$i Harm Phi $k − The fitted value of φ i,k from equation 54. This is output only if one of the keywords "outampphase" or "outampradphase" is given. If "outampphase" is given, then φ i,k /(2π) is reported.
Killharm Per$i Harm R $k 1 − The fitted value of the relative amplitude R i,k,1 from equation 52. This is output only if the "outRphi" or "outRradphi" keyword is given. Killharm Per$i Harm Phi $k 1 − The fitted value of the relative phase φ i,k,1 from equation 56. This is output only if the "outRphi" or "outRradphi" keyword is given. If "outRphi" is given, then φ i,k,1 /(2π) is reported. Killharm Per$i Amplitude − The peak-to-peak amplitude of the full Fourier series (fundamental, sub-harmonics and harmonics) for period $i. Table B.75: Output Columns for "-linfit" Command (Section 2.2.3.3) Linfit $paramname − The fitted value for parameter $paramname, where $paramname is the name of a variable in 'paramlist' supplied on the command line. The Linfit $paramname and Linfit err$paramname columns for each parameter $paramname are output together, followed by the respective columns for the next value of $paramname, until the results have been shown for all varied parameters. Linfit err$paramname − The formal uncertainty for parameter $paramname, this is accurate only if the noise in the light curve is Gaussian and uncorrelated in time, and if the measurement uncertainties supplied in the light curve are correct. Table B.76: Output Columns for "-LS" Command (Section 2.2.1.1) LS Period $peaknum − The period for peak number $peaknum. $peaknum runs from 1 for the highest peak in the periodogram to 'Npeaks', the number of peaks to identify as given on the command line. Note that all columns for a given peak are output together, and then the same quantities are repeated for the next peak, until all 'Npeaks' peaks have been reported. Log10 LS Prob $peaknum − The base-10 logarithm of the false alarm probability for peak number $peaknum, calculated using equation 3, or following Press et al. (1992) if the "noGLS" keyword is given. LS Periodogram Value $peaknum − The value of the Generalized L-S periodogram (Eq. 1) at peak number $peaknum. If the "noGLS" keyword is given, this is instead the traditional normalized L-S periodogram value. LS SNR $peaknum − The spectroscopic S/N ratio for peak number $peaknum (Eq. 6). LS PeriodFix − The fixed period at which the Generalized L-S periodogram statistics are to be evaluated. This is output only if the "fixperiodSNR" keyword is given. Log10 LS Prob PeriodFix − The base-10 logarithm of the false alarm probability at the fixed period. This is output only if the "fixperiodSNR" keyword is given. LS Periodogram Value PeriodFix − The value of the Generalized L-S periodogram at the fixed period. This is output only if the "fixperiodSNR" keyword is given. LS SNR PeriodFix − The spectroscopic S/N at the fixed period. This is output only if the "fixperi-odSNR" keyword is given. − The parameter f 0 in equation 64. This is either the input value if the parameter is not varied (the "novary" keyword was used for this parameter), or the fitted parameter if it is. This also holds for the other parameters output by this command.

Microlens f1
−  This is output only if the "mcmc" fitting mode is used. Here $expression[$i] is the $i-th expression listed in 'exprlist' on the command-line if the "chainstats" keyword is used, or the $ith parameter name in 'paramlist' if the "chainstats" keyword is not used. The term $stats[$j] is a string based on the $j-th statistic in 'statslist' (if the "chainstats" keyword is not given then 'median,stddev' is assumed for 'statslist'). A separate column is output for each combination of $i and $j, with the loop over $j nested within the loop over $i. The possible values for $stats[$j] are: "MEAN", "WEIGHTEDMEAN", "MEDIAN", "WEIGHTED-MEDIAN", "STDDEV", "MEDDEV", "MEDMEDDEV", "MAD", "KURTO-SIS", "SKEWNESS", "PCT$pctval" (where $pctval is the percentile value to two decimal places), "WPCT$pctval", "MAX", "MIN" and "SUM". These terms are output only if the "mcmc" fitting mode is used. The factor by which the magnitude uncertainties are scaled to achieve χ 2 /dof = 1 for that light curve. The minimum JD used for restricting the light curve times. This is output only if the "JDrange" or "JDrangebylc" keyword is used. RestrictTimes MaxJD− The maximum JD used for restricting the light curve times. This is output only if the "JDrange" or "JDrangebylc" keyword is used.   The period for the harmonic series fit simultaneously to the light curve with the transit model (this is included only if the 'fit harm' flag is set to '1'). SoftenedTransit Subharm $k Sincoeff − Coefficient of the sine function for sub-harmonic $k where $k varies from 2 to 'nsubharm'+1 (i.e., this is c 1,$k in Eq. 49 with N P = 1). This is output only if the 'fit harm' flag is set to '1'. SoftenedTransit Subharm $k Coscoeff − Coefficient of the cosine function for sub-harmonic $k where $k varies from 2 to 'nsubharm'+1 (i.e., this is d 1,$k in Eq. 49 with N P = 1). This is output only if the 'fit harm' flag is set to '1'. SoftenedTransit Fundamental Sincoeff − Coefficient of the sine function for the fundamental mode in the harmonic series fit simultaneously to the transit model (i.e., this is a 1,1 in Eq. 49 with N P = 1). This is output only if the 'fit harm' flag is set to '1'.

SoftenedTransit Fundamental Coscoeff −
Coefficient of the cosine function for the fundamental mode in the harmonic series fit simultaneously to the transit model (i.e., this is b 1,1 in Eq. 49 with N P = 1). This is output only if the 'fit harm' flag is set to '1'. SoftenedTransit Harm $k Sincoeff − Coefficient of the sine function for harmonic $k where $k varies from 2 to 'nharm'+1 (i.e., this is a 1,$k in Eq. 49 with N P = 1). This is output only if the 'fit harm' flag is set to '1'. SoftenedTransit Harm $k Coscoeff − Coefficient of the cosine function for harmonic $k where $k varies from 2 to 'nharm'+1 (i.e., this is b 1,$k in Eq. 49 with N P = 1). This is output only if the 'fit harm' flag is set to '1'. SoftenedTransit chi2perdof − χ 2 per degree of freedom for the best-fit model.  A separate column is output for each combination of $i and $j, with the loop over $j nested within the loop over $i.