The VARTOOLS Light Curve Analysis Program

The VARTOOLS program is a command line utility that provides tools for calculating variability/periodicity statistics of light curves as well as tools for modifying light curves. The program is run by issuing a sequence of commands to perform actions on light curves, each command is executed in turn with the resulting light curves passed to the next command. Statistics computed by each command are sent to stdout as an ascii table.

This program is available without any warranty - use it at your own risk. If you do find a bug, please let me know so that I can fix it (jhartman AT cfa DOT harvard DOT edu). Also feel free to let me know if there are any commands/algorithms that you would like to see added to this program.

Citation: If you use this program, please cite Hartman et al., 2008, ApJ, 675, 1254. Also, be sure to cite the relavent sources for the algorithms that you use.

Download Version 1.152 (Aug 13 2009)

Download Version 1.151 (May 19 2009)

Download Version 1.15 (Mar. 3. 2009)

Download Version 1.1 (Dec. 11. 2008)

Download Version 1.0 (Sep 25. 2008)

To install, unpack the program to its own directory. Then type "make" from that directory. Type "./vartools" for the syntax, "./vartools -listcommands" for a list of available commands, or "./vartools -help" for a detailed description. The program is written in ansi c and only uses standard libraries, so it should compile on any linux/unix system without difficulty (let me know though if this is not the case for you). Note that it has been reported to compile successfully on MAC OS X 10.5.6 after deleting the "#include <malloc.h>" statements from the files "commands.h", "detrendtfa.c" and "svdcmp.c" (these includes have been removed in version 1.15).

Commands: In the current version, the following commands are supported:

Calculate the alarm variability statistic for each light curve. Cite Tamuz, Mazeh, and North 2006, MNRAS, 367, 1521 if you use this tool.

Perform an AoV period search on the light curves using phase binning. Specify "Nbin" and a number to change the number of bins used from the default value of 8. minp and maxp are the minimum and maximum periods to search. The intial search will use a frequency resolution of subsample/T where T is the baseline of the light cruve. The peak periods will be refined using a resolution of finetune/T. The program will output the Npeaks highest peaks in the periodogram of the light curve. By default the program will output, for each peak, the period, the theta_aov statistic, the signal to noise ratio (theta_aov - <theta_aov>)/RMS(theta_aov), and the negative natural logarithm of the formal false alarm probability (this is calculated from the value of theta_aov using the Horne and Baliunas 1986, ApJ, 302, 757 estimate for the bandwidth penality). operiodogram should be either 0 or 1. If it is set to 1 then the periodogram for each light curve will also be output to the directory outdir, with the suffix ".aov". The first column in the output is the period, the second column is theta_aov. If the "whiten" keyword is given, then the light curve will be whitened at each peak period and the periodogram will be recomputed before searching for the next peak period. The average and RMS theta_aov used for each peak are computed on the whitened periodogram. The output spectrum will contain the Npeaks computed periodograms. Use the "clip" keyword to change the clipping parameters for calculating the average and RMS of the power spectrum when computing the SNR value of a peak. clip is the sigma-clipping factor, and clipiter is a flag that is 1 or 0 to toggle iterative clipping. By default iterative 5-sigma clipping is used. Use the keyword "uselog" to output the log(theta_aov) SNR instead of the standard statistics (this is the original behavior for this command), i.e. the statistic is (<-ln(theta_aov)> - ln(theta_aov))/RMS(-ln(theta_aov)) where the average and RMS are taken over the entire periodogram for a light curve, in this case the output will also include the average and RMS of -ln(theta_aov). Use the "fixperiodSNR" option to output the AoV statistic, SNR etc at a specified period. See the help for the "-LS" command for an explanation of the syntax. Cite Schwarzenberg-Czerny, A., 1989, MNRAS, 241, 153 and Devor, J., 2005, ApJ, 628, 411 if you use this tool. (The Devor citation is needed because this code is based on his implementation of AoV).

Perform an AoV period search on the light curves using a multi-harmonic model. Nharm is the number of harmonics to use. minp and maxp are the minimum and maximum periods to search. The initial search will use a frequency resolution of subsample/T where T is the baseline of the light cruve. The peak periods will be refined using a resolution of finetune/T. The program will output the Npeaks highest peaks in the periodogram of the light curve. The output for each peak includes the period, the theta_AoV statistic, the signal to noise ratio (theta_aov - <theta_aov>)/RMS(theta_aov), the false alarm probability (see the help for the "-aov command"), and the average and RMS of theta_aov. operiodogram should be either 0 or 1. If it is set to 1 then the periodogram for each light curve will also be output to the directory outdir, with the suffix ".aov_harm". The first column in the output is the period, the second column is theta_aov. If the "whiten" keyword is given, then the light curve will be whitened at each peak period and the periodogram will be recomputed before searching for the next peak period. The average and RMS theta_aov output for each peak are computed on the whitened periodogram. The output spectrum will contain the Npeaks computed periodograms. Use the "clip" keyword to change the clipping parameters for calculating the average and RMS of the power spectrum when computing the SNR value of a peak. clip is the sigma-clipping factor, and clipiter is a flag that is 1 or 0 to toggle iterative clipping. By default iterative 5-sigma clipping is used. Use the "fixperiodSNR" option to output the AoV statistic, SNR etc at a specified period. See the help for the "-LS" command for an explanation of the syntax. Cite Schwarzenberg-Czerny, A., 1996, ApJ, 460, L107 if you use this tool.

Calculate the discrete auto-correlation function (Edelson and Krolik 1988, ApJ, 333, 646), these are written out to outdir/basename.autocorr, start, stop and step are the times in days for sampling the auto-correlation. Note that rather than using the variance of the light curve in the denominator, we use the formal uncertainty (we do not subtract the "measurement error" from the variance as done in the Edelson and Krolik formula since this could lead to imaginary numbers in the case where measurement errors are overestimated). If you wish to use the variance in the denominator rather than the formal uncertainty you should issue the -changeerror command before calling this routine. Note that due to binning, when the variance is used in the denominator the autocorrelation function may be smaller than 1 unless the time step used is less than the time difference between any consecutive measurements in the light curve.

Bin the light curves in time (or phase if a -Phase command has already been issued). combinetype is 0 to take the average of points in a bin, 1 to take the median of the points or 2 to take the weighted average. One should specify either the binsize in units of the time coordinate (e.g. in phase if the light curves have been phased), or the number of bins to split the light curve time-span into. By default the first bin begins at the initial time in the light curve (t0), if firstbinshift is specified, then this will be shifted by t0 - firstbinshift/binsize. timetype should be 0, 1 or 2. If it is 0 then the output time for each bin is the time at the center of the bin, if it is 1 then the time will be the average of the times of points that fall within the bin, if it is 2 then it will be the median.

This command runs BLS on the light curves. You may either set a fixed minimum and maximum q (fraction of orbit in transit) for the search, or you can specify a minimum and maximum stellar radius to consider (in solar radii) in which case the qmin and qmax values are calculated for each trial period P as q = 0.076 * R**(2/3) / P**(2/3) (with P in days - this assumes R = M as an approximation for the lower main sequence). The minimum and maximum periods to search are given in days, nfreq is the number of trial frequencies (10000 is a decent number). nbins is the number of phase bins to break the light curve into (200 is a decent number). timezone is the number to add to the julian date to get the local date (-7 for arizona), this is used to determine which nights different observations come from and thus to determine the fraction of delta-chi2 that comes from a single night. Npeak is the number of peaks in the BLS spectrum to find and report. outperiodogram is a flag that is set to 1 to output the BLS period vs. SN spectrum. outdir is the output directory for the BLS spectrum if outperiodogram is set to 1, ".bls" will be appended to the filename. omodel is a flag set to 1 or zero that can be used to output the model for the light curve, the output directory is then given in modeloutdir, the suffix ".bls.model" will be appended to the filename. correctlc is either 0 or 1, set to 1 it will subtract the transit model from the light curve before passing it to the next command. Cite Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369 if you use this tool.

This command runs BLS at a fixed period on the light curves - it searches merely for the most transit-like signal at the specified period. The period comes either from the last aov command, from the last ls command, or is specified (in which case it must be in the appropriate column of the lc_list). Like BLS you can either specify a minimum and maximum stellar radius to consider, or a minimum and maximum q to consider. Similarly you must specify the number of bins to use and the timezone. omodel is 1 to output the model light curve to model_outdir (with ".blsfixper.model" appended to the end of the file name. correctlc is 1 to subtract the transit model from the light curve before passing it to the next command. Cite Kovács, G., Zucker, S., & Mazeh, T. 2002, A&A, 391, 369 if you use this tool.

Replace the formal errors in a light curve with the RMS of the light curve.

Calculate chi2 per dof for the light curves. The output will include chi2 and the error weighted mean magnitude.

Calculate chi2 per dof after applying a moving mean filter to the light curves. (Note that the light curves passed to the next command are unchanged). Nbin filters are used (with Nbin resulting estimates of chi2 and the error weighted mean). The width of each filter is given by 2.0*bintime (bintime is given in minutes)

Sigma-clip the light curves using a clipping factor of sigclip. iter is a flag that is either 1 for iterative clipping, or 0 for single pass clipping. The output table will include the number of points that were clipped.

Decorrelate the light curves against specified signals. The signals are either global signals (e.g. the airmass) that are input in separate files (with the format: JD signal_value), or are light curve specific signals (e.g. the sub-pixel position) that are input as additional columns in the light curves. correctlc is a flag that is 0 or 1, if the value is 1 then the light curves passed to the next command will be decorrelated, if the value is 0 then resulting chi2 from decorrelation and the decorrelation coefficients will be output to the table, but the light curves themselves will not be decorrelated. zeropointterm is a flag that is 0 or 1, if the value is 1 then a zeropoint term is included in the decorrelation, if it is 0 then no such term is included. subtractfirstterm is a flag that is 0 or 1, if the value is 1 then the light curves are decorrelated against the signal minus the first signal value (this is useful if one is decorrelating against the JD, for example, in which case one would decorrelate against JD - JD0 where JD0 is the first JD in the light curve), if it is 0 then the light curves are decorrelated against the signal without subtracting the first term. Nglobalterms is the number of global files. globalfile1...globalfileN or the names of those files, order1...orderN are the orders of the polynomials used to decorrelate each signal (each must be greater than 0). Nlcterms is the number of light curve specific signals. The columns of these signals are given by lccolumn1...lccolumnN. The orders of the polynomials are given by lcorder1...lcorderN. omodel is a flag set to 1 or zero that can be used to output the model for the light curve, the output directory is then given in modeloutdir, the suffix ".decorr.model" will be appended to the filename.

Compute the DFT power spectrum of the light curve and optionally apply the clean algorithm (Roberts et al. 1987) to it. nbeam is the number of points per 1/T frequency element to include in the power spectrum. Use "maxfreq" to set the maximum frequency to maxf. The default value for the maximum frequency is 1/(2*min_time_separation). Use "outdspec" to output the dirty power spectrum, the suffix ".dftclean.dspec" will be appended to the name of the output file. Use "finddirtypeaks" to find the top Npeaks peaks in the dirty power spectrum. Use "outwfunc" to output the window function, with the suffix ".dftclean.wfunc". Use "clean" to apply the clean algorithm to the power spectrum. Please cite Roberts, D.H., Lehar, J., and Dreher, J.W. 1987, AJ, 93, 4 if you use this algorithm. To use clean you must specify the gain which is a value between 0.1 and 1 (with lower values giving slower convergence). The procedure will continue to clean the spectrum until the last peak is less than SNlimit times the noise. Use "outcbeam" to output the clean beam which is convolved with the deconvolved spectrum to produce the final clean power spectrum. The suffix for the clean beam file is ".dftclean.cbeam". Use "outcspec" to output the clean power spectrum, the suffix is ".dftclean.cspec". Finally use "findcleanpeaks" to find the top Npeaks peaks in the clean power spectrum. This routine uses the FDFT algorithm by Kurtz 1985, MNRAS, 213, 773 to compute the DFT.

Rescale the magnitude uncertainties by fitting a linear relation between (expected rms)^2 vs (chi2)*(expected rms)^2 for all light curves so that chi2 per dof is distributed about unity, Outliers in the distribution are clipped using a clipping factor of sigclip. This routine requires a list of light curves. If this command is invoked, then all light curves will be read into memory. The output will include the average rescale factor for each light curve (taken to be sqrt(chi2_after/chi2_before)).

Use this routine to find variability blends. The routine operates by matching a list of potential variables to a list of light curves. For each potential variable, the flux amplitude of all matching light curves will be measured and the one with the highest amplitude will be chosen. The name and flux amplitude of the highest amplitude light curve are included in the output statistics table. If this routine is used, light curves must be read from an input list with x and y coordinates given as additional columns in the list file (i.e. you must use "-l" for input rather than "-i"), this will be used as the list of potential variables. matchrad is the matching radius used to match stars that are possibly blended. If "radec" is specified, then matchrad is given in arcseconds and the x and y coordinates are ra/dec in degrees, if "radec" is not specified then rectangular matching will be used. Use the "fix", "spec" or "fixcolumn" keywords to specify the source for the periods to use for each potential variable. If "fix" is given then the periods for all potential variables are set to the listed value. If "spec" is given then the periods should be included as an additional column in the input list file (after the x and y coordinate columns). If "fixcolumn" is given then the period for each potential variable in the list will be set to an arbitrary previously computed statistic by giving the name or number of the output column as seen with the -header or -numbercolumns options. By default the list of potential variables will be matched to itself. If, however, "starlist" is specified, then the list of potential variables will be matched to the list contained in the file starlistfile. The first three columns in this file are the light curve name, and the x and y coordinates. Use "zeromag" to specify the zero-point magnitude for converting magnitudes into fluxes, the default value is 25.0. If you do not wish to convert from magnitudes into fluxes (e.g. if the input light curves are already in flux units), use "nofluxconvert", make sure you know what you're doing then when you use this option. The amplitude of the light curves are taken to be the peak-to-peak amplitude of a fourier series fit to the light curve with period P. If the period is less than or equal to zero then the amplitude will be measured by fitting the fourier series to the unphased light curve (i.e. the period will be set equal to the total duration of the light curve). Use "Nharm" to specify the number of harmonics to include in the Fourier series (if it is 0 then only a sinusoid will be fit, the default value is 2). You can output the names and flux amplitudes of all stars matching to each potential variable by using "omatches" and giving the name of the file to output this information to.

Convert light curves from ISIS differential flux into magnitudes. In this case a light curve list must be used, and the reference magnitudes of the stars (after aperture correction) should be given as an additional column in the light curve list. mag_constant is the magnitude of a source with a flux of 1 ADU, offset is an additive constant to apply to the output light curves. This command does not yield any output to stdout, you must explicitly call -rms or -chi2 if you want statistics.

This routine determines the minimum amplitude that a light curve could have and still be detected at a given period (which is either specified in the input list, or taken from the last LS command) with a LS log(formal false alarm probability) < thresh. The light curve signal is calculated by fitting a fourier series at the period of the light curve with Nharm harmonics and Nsubharm subharmonics if "harm" is specified, or it is read in from a file if "file" is specified. If reading the signal from a file, listfile is a file with two columns of the form: signal_file signal_amp, with one line for each light curve being operated on. Each signal file should contain the signal magnitude in the third column, signal_amp is the amplitude in magnitudes of the signal (this is to allow for cases where the signal amplitude is greater than just the max observed signal magnitude minus the min observed signal magnitude, e.g. if the signal comes from fitting a fourier series to the light curve at a different period from the one used to do the ls selection). To calculate the minimum amplitude, the signal is subtracted it from the light curve and re-added after scaling it by a number. minP is the minimum period that would be searched (this is needed since it sets the false alarm probability scale). The output is the minimum factor by which the signal could be scaled and still be detectable together with the minimum peak-to-peak amplitude.

inject a harmonic series into the light curves. The harmonic series has the form: A_1*cos(2*pi*(t/P + phi_1)) + sum_{k=2}^{Nharm+1} A_k*cos(2*pi*(t*k/P + phi_k)) + sum_{k=2}^{Nsubharm+1} A_k*cos(2*pi*(t/k/P + phi_k)). The period for the series is either specified for each light curve in the list file (use "spec"), it is fixed to a particular value for all light curves (use "fix"), it is generated randomly either from a uniform distribution (use "rand") or from a uniform log distribution (use "logrand"), or the frequency is generated from a uniform distribution (use "randfreq") or from a uniform log distribution (use "lograndfreq"). Specify the number of harmonics, and then for each harmonic specify how the amplitude and phase of that harmonic is to be generated. Note that the fundamental period in this case corresponds to harmonic 1, so for Nharm=0 there must be one set of amplitude/phase commands, for Nharm=1 there are two sets, etc. ampspec, ampfix, amprand and amplogrand specify whether the amplitude for the harmonic is specified as a column in the light curve list file, or if it is fixed to a particular value for all light curves, or if it is randomly generated from a uniform or uniform log distribution. You may optionally specify "amprel" after the amplitude command to make the specified amplitude of this harmonic to be relative to the amplitude of the fundamental mode (i.e. the input is R_i1 = A_i/A_1 rather than A_i). This would be used, for example, if one wishes to simulate a specific signal form (say a Cepheid) but with a random over-all amplitude. phasespec, phasefix and phaserand specify whether the phase of the harmonic is specified as a column in the light curve list file, or if it is fixed to a particular value for all light curves, or if it randomly generated from a uniform distribution between 0 and 1. The phase is the phase at time t=0. If the optional "phaserel" term is given then the phase for the harmonic will be relative to the phase of the fundamental (i.e. the input is phi_k1 = phi_k - k*phi_1 rather than phi_k). Then specify the number of sub-harmonics and list how the amplitudes and phases for each sub-harmonic are to be generated (periods k*P with k=2...Nsubharm+1). omodel is a flag that is either 1 or 0 depending on whether or not the model light curve will be output, if it is 1 then the output directory is given in modeloutdir, the suffix ".injectharm.model" will be appended to the filename.

Inject a Mandel-Agol transit model into the light curve. The source for the following parameters must be given: period in days (or freq in 1/day), radius of the planet in jupiter radii, mass of the planet in jupiter masses, phase of the transit at time T=0 (phase = 0 corresponds to transit center), sin_i of the transit, e/omega (omega in degrees) or h/k=esinomega/ecosomega, mass of the star in solar masses, radius of the star in solar radii and the limb darkening coefficients. Each parameter can either be specified as a column in the input list (use the *spec keywords), or fixed to a specific value given on the command line (use the *fix keywords). In some cases the parameter can be generated from a uniform random or uniform log random distribution. If the "sinirand" keyword is used then sin_i is drawn from a uniform orientation distribution with the constraint that there must be a transit (for a circular orbit this corresponds to cos(i) being uniform between 0 and (Rstar + Rp)/a). omodel is a flag that is either 1 or 0 depending on whether or not the model light curve will be output, if it is 1 then the output directory is given in modeloutdir, the suffix ".injecttransit.model" will be appended to the filename.

Calculate Stetson's J statistic, L statistic and the Kurtosis of each light curve. The timescale is the time in minutes that distinguishes between "near" and "far" observations. The dates file should contain JDs for all possible observations in the first columns - this is used to calculate the maximum possible weight. Cite Stetson, P.B. 1996, PASP, 108, 851 if you use this tool.

This command whitens light curves against one or more periods. The mean value of the light curve, the period of the light curve and the cos and sin coefficients are output. The light curves passed to the next command are whitened light curves (unless the keyword "fitonly" is given). The origin of the period(s) is either from the most recent previous aov command (either aov, or aov_harm), from the most recent previous LS command, two periods one from aov, the other from LS, the period from the most recent injectharm command, Nper periods per1 through perN that are fixed for all light curves, or Nper periods specified in the input list (the periods are read off in order as additional columns in the input light curve list - a list must be used for this option). The light curves will be whitened using Nharm higher-harmonics (frequencies of 2*f0, 3*f0, ... (Nharm + 1)*f0) and Nsubharm sub-harmonics (frequencies of f0/2, f0/3, ... f0/(Nsubharm + 1)). omodel is a flag set to 1 or zero that can be used to output the model for the light curve, the output directory is then given in modeloutdir, the suffix ".killharm.model" will be appended to the filename. If "fitonly" is specified, then the model is not subtracted from the light curve (a keyword is used rather than a flag to maintain compatability with scripts written before this option was added). By default the a_k and b_k cos and sin coefficients are output. If the keyword "outampphase" or "outampradphase" is given, then the output will be the amplitudes A_k = sqrt(a_k^2 + b_k^2) and the phases (phi_k = atan2(-b_k, a_k)/2pi for "outampphase" or phi_k = atan2(-b_k, a_k) for "outampradphase"). If the keyword "outRphi" or "outRradphi" is given then the output will be the relative amplitudes R_k1 = A_k/A_1 and phases phi_k1 = phi_k - k*phi_1 (in units of 0 to 1, or in radians for the two keywords respectively). Note that for sub-harmonics, k = 1/2, 1/3, etc. For the fundamental mode the amplitude A_1 and the phase phi_1 will be given.

Perform a Lomb-Scargle period search on the light curves. In this case the program outputs the signal to noise ratio (SNR) and the formal false alarm probability for the peaks (assuming white noise with rms equal to the rms of the light curve). The output periodogram has the format: frequency normalized-LS statistic. This version does not allow the mean to float. If the "whiten" keyword is given, then the light curve will be whitened at each peak period and the periodogram will be recomputed before searching for the next peak period. If "whiten" is used then the RMS used in finding the signal-to-noise ratio is computed on the whitened periodogram. The output spectrum will contain the Npeaks computed periodograms. Use the "clip" keyword to change the clipping parameters for calculating the average and RMS of the power spectrum when computing the SNR value of a peak. clip is the sigma-clipping factor, and clipiter is a flag that is 1 or 0 to toggle iterative clipping. By default iterative 5-sigma clipping is used. If the "fixperiodSNR" keyword is given, then the FAP and SNR will be given for a specific period as well as for the peaks. This period can either be from the most recent issued aov or aov_harm command, from the most recent ls command, from the most recent injectharm command, it can be fixed to a particular value for all light curves, it can be specified as an additional column in the input list, or it can be set to an arbitrary previously computed statistic (using the "fixcolumn" option together with the name or number of the output column as seen with the -header or -numbercolumns options). Cite Lomb, N.R. 1976, A&SS, 39, 447, Scargle, J.D. 1982, ApJ, 263, 835, Press, W.H. & Rybicki, G.B. 1989, ApJ, 338, 277, and Press, W.H., Teukolsky, S.A., Vetterling, W.T. & Flannery, B.P. 1992, Numerical Recipes in C, 2nd ed. (New York: Cambridge University Press) if you use this tool.

Fit a Mandel and Agol (2002) transit model to the light curve. The initial parameters either are automatically estimated using the results from bls if "bls" is specified, a blsfixper command, or they are entered directly in the command line as the values for P0 T00 etc. The parameters are P0 - period, T00 - time of transit, r0 - ratio of planet radius to star radius, a0 - ratio of semi-major axis to star radius, sin_i0 - sine of the inclination angle, e0 - eccentricity, omega0 - argument of periastron in degrees, and mconst0 - the out of transit magnitude. If mconst0 is negative then it will be estimated directly from the light curve. "quad" or "nonlin" are used to specify the limb darkening law, either quadratic or non-linear, following Claret. The number of initial limb darkening coefficients must be 2 for quadratic limb-darkening or 4 for non-linear limb-darkening. fitP fitT0 ... should be 1 or 0 depending on whether or not the corresponding parameter is allowed to vary. If fitRV is 1 then the program will simultaneously fit an RV curve stored in the file RVinputfile (first column JD, second column RV, third column RV error). The model RV will be output to the file RVmodeloutfile (this will be a model evaluated at a number of evenly spaced phase points rather than at the observed RV phases). If fitting RV specify the initial K0 and gamma0 together with flags fitK and fitgamma to specify whether or not these parameters are allowed to vary in the fit. Set correctlc to 1 to subtract the model from the light curve, set it to 0 to only perform the fit. omodel is a flag set to 1 or zero that can be used to output the model for the light curve, the output directory is then given in modeloutdir, the suffix ".mandelagoltransit.model" will be appended to the filename. Cite Mandel, K., & Agol, E. 2002, ApJ, 580, L171 if you use this tool.

Fit a simple microlensing model to the light curve. This program uses the functional form for the model given by Wozniak, P. R., et al. 2001, AcA, 51, 175. For each parameter you can optionally specify the source for the initial parameter value: either "fix" the initial value for all light curves to the value fixval, use "spec" to specify the value as an additional column in the input list, use "fixcolumn" to take the initial value from a previously computed statistic, or use "auto" to automatically determine the initial parameter value. The default for each parameter is "auto". You can then optionally specify an initial step size for each parameter for use in the downhill simplex fit. You may also optionally use "novary" to not vary this parameter in the fit. Use the "correctlc" keyword to subtract the model fit from the light curve before passing to the next command, and the "omodel" keyword to output the model. The suffix will be "microlens".

output the light curves to directory outdir

Replace the time variable of a light curve with its phase and sort the light curve by the phase. The period for phasing is either taken from the last aov command, the last ls command, the last BLS command, or is specified in the list (in this case the list must be used).

Rescale the magnitude uncertainties of each light curve so that Chi2 per dof is equal to 1 for every light curve. The rescale factor for each light curve will be included in the output table.

This command can be used in conjunction with the -savelc command to restore the light curve to a previous state. savenumber is 1, 2, 3, ... to restore the light curve to its state at the first, second, third ... call to -savelc. For example, suppose you want to try running TFA followed by LS on the light curve using 3 different template lists. A command string of the form: "... -savelc -TFA trendlist1 ... -LS ... -restorelc 1 -TFA trendlist2 ... -LS ... -restorelc 1 -TFA trendlist3 ... -LS ..." would accomplish this as each time "-restorelc 1" is called the light curve is restored to its form at the first -savelc command.

Calculate the rms of the light curves. The output will include RMS, the mean magnitude, the expected RMS and the number of points in the light curve

Similar to chi2bin, this calculates RMS after applying a moving mean filter.

This command can be used in conjunction with the -restorelc command to save a light curve and later restore it to this state. For example, suppose you want to try running TFA followed by LS on the light curve using 3 different template lists. A command string of the form: "... -savelc -TFA trendlist1 ... -LS ... -restorelc 1 -TFA trendlist2 ... -LS ... -restorelc 1 -TFA trendlist3 ... -LS ..." would accomplish this as each time "-restorelc 1" is called the light curve is restored to its form at the first -savelc command.

Fit a Protopapas, Jimenez and Alcock transit model to the light curve. The initial parameters either come from bls if "bls" is specified, a blsfixper command if "blsfixper" is specified, otherwise they are entered directly in the command line as the values for P0 T0 etc. If mconst0 is < 0 then it will be estimated directly from the light curve. fitP fitT0 should be 1 or 0 depending on whether or not the corresponding parameter is allowed to vary. Set correctlc to 1 to subtract the model from the light curve, set it to 0 to only perform the fit. omodel is a flag set to 1 or zero that can be used to output the model for the light curve, the output directory is then given in modeloutdir, the suffix ".softenedtransit.model" will be appended to the filename. If fit_harm is 1 then simultaneously fit a harmonic series to the light curve. The period comes from either the last aov, the last ls, the last bls model or is specified in the command line as Pharm (you must type "spec" and then the value), finally the number of harmonics and sub-harmonics must also be specified. Cite Protopapas, P., Jimenez, R., & Alcock, C. 2005, MNRAS, 362, 460 if you use this tool.

This command fits a single, circular, uniform temperature spot model to the light curve using the Dorren 1987 model. Use aov or ls to choose the initial period to be from AoV or LS, use spec to specify the initial period in the input list (the period is read off as an additional column in the input light curve list - a list must be used for this option). The other parameters are initial guesses as in Dorren 1987, mconst is the constant magnitude term. Set mconst0 to a negative value to have the program automatically estimate the initial mconst. Note that here we take a = a_d/(pi * (1 - mu_* / 3)) and b = b_d/(pi * (1 - mu_* / 3)) where a_d and b_d are the a and b terms from Dorren 1987. fitP... are flags denoting whether or not a parameter is to be varied. Set correctlc to 1 to subtract the model from the light curve, set it to 0 to only perform the fit. omodel is a flag set to 1 or zero that can be used to output the model for the light curve, the output directory is then given in modeloutdir, the suffix ".starspot.model" will be appended to the filename. Cite Dorren 1987, ApJ, 320, 756 if you use this tool.

Run the SYSREM algorithm to identify and remove trends from an ensemble of light curves. If this command is used, the -readall option will automatically be set. Also, if this command is used then one must have as input a light curve list. A total of Ninput_color + Ninput_airmass trends will be searched for. Ninput_color of these will have the colors set initially, Ninput_airmass will have the airmass terms set initially. The initial color terms will be read in as additional columns in the input light curve list, there must be Ninput_color of these columns. The initial_airmass_file gives a list of the initial airmass trends to use, the first column is JD and the subsequent Ninput_airmass columns are the initial airmass trends. As implemented this routine first filters the trends with the airmass terms specified initially and then filters the trends with the color terms specified initially. sigma_clip1 is the sigma clipping used in calculating the mean magnitudes for the light curves. sigma_clip2 is the sigma clipping used in determining whether or not points contribute to the airmass or color terms when doing the fit. Any points with magnitude less than saturation will not contribute. omodel is a flag that is 1 or zero for outputing the model light curves, the output directory is given in model_outdir, the suffix ".sysrem.model" will be appended to the filename. The format for the model light curve is: JD mag mag_model sig clip, where clip is either 1 if the point was included in determining the trends or 0 if it wasn't. otrends is a flag that is 1 or zero for outputing the final trends. These will be output to the file trend_outfile, the first column is the JD and subsequent columns are for each trend signal. If you use this command, be sure to cite Tamuz, Mazeh and Zucker (2005, MNRAS, 356, 1466).

Run the Trend Filtering Algorithm on the light curves. The trendlist is a file containing a list of light curves to be used as trends, it has the format: trendname trendx trendy, where trendname is the file name and trendx and trendy are the coordinates of the trend star. You can optionally specify the readformat for the light curves in the trendlist, by default Nskip is 0, jdcol is 1 and magcol is 2. Note that if -matchstringid is set then the JDs of the trend light curves are not read-in and instead jdcol corresponds to the column storing the string-ids for point matching between light curves. dates_file is a file giving all the dates in the second column, the unique string-ids are read-in from the first column if the -matchstringid option is used (it is the same format as the dates file for the ISIS image subtraction program where the first column is the filename and the second is JD). Trend stars within pixelsep of the light curve in question will not be used in detrending the light curve. To you this routine you need an input light curve list, the x and y positions of each light curve must be given as additional columns in the list. correctlc is a flag that is 0 or 1 indicating whether the light curves passed to the next command should have the tfa applied. ocoeff is a flag to denote whether or not to output the list of coefficients for each trend for a given light curve, if set to 1 then they will be output to coeff_outdir with ".tfa.coeff" appended to the end of the filename. omodel is a flag set to 1 or zero that can be used to output the tfa model for the light curve, the output directory is then given in model_outdir, the suffix ".tfa.model" will be appended to the filename. If you use this routine be sure to cite Kovacs, Bakos and Noyes, 2005, MNRAS, 356, 557

Run the Trend Filtering Algorithm in signal reconstruction mode on the light curves (i.e. iteratively filter the light curve and fit a simple signal to the light curve). In addition to the parameters used by the -TFA command, this command allows for simulataneously decorrelating the light curve against additional light curve specific signals. This is useful, for example, if one wishes to do EPD (external parameter decorrelation) and TFA on a high-amplitude variable star. To do this, specify "decorr", the iterativeflag is 1 if the decorrelation and the TFA will be done iteratively (this is faster) or 0 if they will be done simultaneously (more correct, but slower). The Nlcterms, lccolumn and lcorder terms are the same as for the -decorr command. Any global signals to decorrelate against should just be included in the TFA trendlist. Other parameters that are different from the -TFA command include: "dotfafirst" is a flag that is 0 or 1, if the flag is 1 then TFA will be applied to the input light curve first and the signal will then be determined on the residual in each iteration, if it is 0 then the signal is determined and subtracted from the light curve and TFA is applied to the residual in each iteration. The iterations will stop once the fractional change in the RMS is less than "tfathresh" or if the number of iterations reaches "maxiter". The model signal is either taken to be the mean value of the binned light curve (specify "bins" and then the number of bins to use), it is a fixed signal form read in from a file (specify "signal" and then a filename), or it is a Fourier series that is simultaneously fit to the light curve with TFA (specify \"harm\"). If binning is used, then you can specify "period" to make the binning phase-binning, and then specify where the period should be read from. If "spec" is used then the period is read from the input light curve list (it should be put after the x and y coordinates for the star), if "fix" is used then the period is fixed to the specified period value for all light curves. If "signal" is used, then the filename will be a file containing a list of signal files, one file for each light curve, with each signal file containing the signal in the second column. The quantity a*signal + b is fit to the light curve, where a and b are free parameters. If "harm" is used, then the signal will be a Fourier series that is fit simultaneously to the light curve with TFA. In this case there is no TFA iteration. Nharm is the number of harmonics in addition to the fundamental to include in the fourier series, Nsubharm is the number of subharmonics. If "period" is specified then the period for the fourier series will be taken from the specified source. If "period" is not specified then the period for the Fourier series will be set to the time-span of the observations. If you use this routine be sure to cite Kovacs, Bakos and Noyes, 2005, MNRAS, 356, 557.

either provide an individual light curve, or a list of light curves. If lcname or lclist is "-" then input will be taken from stdin. Some commands can only be used with lists. The list should contain a single light curve filename per line in the first column. Additional columns may be necessary to specify, for example, the periods to use for pre-whitening or any other values required by some commands.

use this option to provide a one line header at the start of the output table.

use this option to output the header and then quit without processing any light curves.

use this option to use a tab-delimited starbase format for the output table rather than the default space-delimited format.

use this option to print only the basename of each light curve in the output table.

use this option to specify the format of the input light curves. Nskip is the number of lines to skip at the beginning of each file (this includes any lines that begin with '#' which are otherwise automatically ignored), the default value is 0. If you need to read in a column of strings from each light curve specify "stringid" and then the column, by default no column of strings is read in. The stringid column must be specified in this fashion, however, if the matchstringid option is set. col_jd, col_mag and col_sig are the columns that contain the time, magnitude and magnitude uncertainties (or differential flux and differential flux uncertainty if the light curve will be passed through the -fluxtomag command), the default values are 1, 2 and 3.(Note that the time is assumed to be in days).

use this option to force the program to read in all the light curves at once. If not specified, this mode will only be used if a command that requires storing all the light curves in memory is issued.

If this is specified then periodograms will be printed out in ascii format. Otherwise they will be stored in binary.

Time measurements within jdtol of each other are considered equal. The default value is 0.000010 days.

If this option is specified then commands that require image-point matching will use a string-id for each image rather than by comparing the jd values. If this option is set then the "stringid" column for each light curve must be specified with the -readformat option.

If this is specified then stdout will not be buffered, by default it is buffered.

If this is specified the output table will not be generated, however any calls to output light curves, periodograms etc. will still be executed.

Use this to set the seed for the random number generator. The seed should be an integer, or you may use the word "time" to seed the generator with the system clock. If this is not specified the seed value will be 1.

Prefix each column name in the header with its column number.

$ ./vartools -l EXAMPLES/lc_list -rms -decorr 1 1 1 0 1 1 2 0 -rms -tab

Name Mean_Mag_0 RMS_0 Expected_RMS_0 Npoints_0 Decorr_constant_term_1 Decorr_constant_term_err_1 \
LCColumn_1_coeff_1_1 LCColumn_1_coeff_err_1_1 LCColumn_1_coeff_2_1 LCColumn_1_coeff_err_2_1 Decorr_chi2_1 \
Mean_Mag_2 RMS_2 Expected_RMS_2  Npoints_2
---- ---------- ----- -------------- --------- ---------------------- -------------------------- \
-------------------- ------------------------ -------------------- ------------------------ ------------- \
---------- ----- --------------  ---------
EXAMPLES/1 10.24745 0.15944 0.00101 3122 10.0830375984825 0.0000325849746 0.0097933162509 0.0000059117875 \
0.0002554062775 0.0000001956124 6.68601 10.24728 0.00211 0.00101 3122
...

Fit a quadratic polynomial to the light curves given in the file EXAMPLES/lc_list. To do this we use no global terms and 1 lc term. For the lc term we use the first column in the light curve (the JD) and fit a second order polynomial in this term to each light curve. We fit for the zero-point term and correct the light curve (so that commands after -decorr will receive light curves with the best-fit quadratic polynomial removed, note that when the light curve is corrected the mean is kept constant), we also subtract the first term in the signal that we are decorrelating against (use JD - JD0 rather than JD, since JD*JD runs into round-off problems whereas (JD - JD0)*(JD - JD0) does not), but we do not output the corrected light curves to the disk. The rms is determined before and after the fit, the -tab option causes the output table to be in starbase format with headers for each column. To interpret the output, note that for light curve 1 we find that the best-fit quadratic has the form: 0.0002554062775*(JD - 53725.173920)*(JD - 53725.173920) + 0.0097933162509*(JD - 53725.173920) + 10.0830375984825, and that fitting this equation reduces the RMS from 0.15944 mag to 0.00211 mag (a quadratic signal was injected into this particular light curve).

$ ./vartools -i EXAMPLES/2 -LS 1.0 2.0 0.01 1 0 -Killharm ls 0 0 1 EXAMPLES/OUTDIR1 -tab

Name LS_Period_1_0 Log10_LS_Prob_1_0 LS_SNR_1_0 Killharm_Mean_Mag_1\	
Killharm_Period_1_1 Killharm_Per1_Fundamental_Sincoeff_1 Killharm_Per1_Fundamental_Coscoeff_1\
Killharm_Per1_Amplitude_1
---- ------------- ----------------- ---------- -------------------\
------------------- ------------------------------------ ------------------------------------\
-------------------------
EXAMPLES/2 1.23588006 -707.23302 4.68589 10.12178 1.23588006 0.02004 -0.04580 0.09998

Run the Lomb-Scargle period search on the light curve EXAMPLES/2. This light curve has had a sine-curve with a period of 1.2354 days and an amplitude of 0.05 mag injected into it. After using L-S to find the period, fit a harmonic series with no harmonics or sub-harmonics (i.e., just a sine-curve at the period found by L-S) to the light curve. Output the light curve with the best-fit model in the third column to the file EXAMPLES/OUTDIR1/2.killharm.model. Note that in this case L-S finds the period of 1.23588006 days, and the killharm routine fits the function 0.02004*sin((JD - JD0)*2*pi/1.23588006) - 0.0458*cos((JD - JD0)*2*pi/1.23588006) + 10.12178 to the light curve. The killharm routine returns a peak-to-peak amplitude of 0.09998 mag.

$ ./vartools -i EXAMPLES/3 -MandelAgolTransit 2.12345 53725.174 0.1 10.0 1.0 0. 0. 0 quad 0.236 0.391 \
  0 0 0 0 0 0 0 0 0 0 0 0 1 EXAMPLES/OUTDIR1 -tab

Name MandelAgolTransit_Period_0 MandelAgolTransit_T0_0 MandelAgolTransit_r_0 MandelAgolTransit_a_0 \
MandelAgolTransit_sin_i_0 MandelAgolTransit_e_0 MandelAgolTransit_omega_0 MandelAgolTransit_mconst_0 \
MandelAgolTransit_ldcoeff1_0 MandelAgolTransit_ldcoeff2_0 MandelAgolTransit_chi2_0
---- -------------------------- ---------------------- --------------------- --------------------- \
------------------------- --------------------- ------------------------- -------------------------- \
---------------------------- ---------------------------- ------------------------
EXAMPLES/3 2.12345000 53725.17400000 0.10000 10.00000 1.00000 0.00000 0.00000 0.00000 0.23600 \
0.39100 129730808.07162

$ gawk '{print $1, $2 + $3, $4}' EXAMPLES/OUTDIR1/3.mandelagoltransit.model > EXAMPLES/3.transit

These two commands effectively inject a Mandel and Agol transit into a light curve. The first line calls the MandelAgolTransit vartools command to create a simulated transit with a period of 2.12345 days, T0 = 53725.174, Rplanet/Rstar = 0.1, a/Rstar = 10.0, sini = 1.0, a circular orbit, and quadratic limb darkening for the star with coefficients of 0.236 and 0.391. We do not fit for any of the parameters (10 zeros), nor do we fit an RV curve (the 11th 0), nor do we correct the light curve (the 12th 0). The model is output in the third column of EXAMPLES/OUTDIR1/3.mandelagoltransit.model. In the second line we use gawk (or awk) to add the model transit signal to the light curve and output the result to the file EXAMPLES/3.transit. Note that this procedure may now also be done using the -Injecttransit command so that the external call to gawk is not needed.

$ ./vartools -l EXAMPLES/lc_list_tfa -rms -TFA EXAMPLES/trendlist_tfa EXAMPLES/dates_tfa 25.0 1 0 0 \
  -BLS q 0.01 0.1 0.5 5.0 2000 200 7 2 0 1 EXAMPLES/OUTDIR1 1 -rms -tab

Name Mean_Mag_0 RMS_0 Expected_RMS_0 Npoints_0 TFA_MeanMag_1 TFA_RMS_1 BLS_Period_1_2 BLS_SN_1_2 BLS_SR_1_2 \
BLS_SDE_1_2 BLS_Depth_1_2 BLS_Qtran_1_2 BLS_i1_1_2 BLS_i2_1_2 BLS_deltaChi2_1_2 BLS_fraconenight_1_2 \
BLS_Npointsintransit_1_2 BLS_Ntransits_1_2 BLS_Npointsbeforetransit_1_2 BLS_Npointsaftertransit_1_2 \
BLS_Rednoise_1_2 BLS_Whitenoise_1_2 BLS_SignaltoPinknoise_1_2 BLS_Period_2_2 BLS_SN_2_2 BLS_SR_2_2 \
BLS_SDE_2_2 BLS_Depth_2_2 BLS_Qtran_2_2 BLS_i1_2_2 BLS_i2_2_2 BLS_deltaChi2_2_2 BLS_fraconenight_2_2 \
BLS_Npointsintransit_2_2 BLS_Ntransits_2_2 BLS_Npointsbeforetransit_2_2 BLS_Npointsaftertransit_2_2 \
BLS_Rednoise_2_2 BLS_Whitenoise_2_2 BLS_SignaltoPinknoise_2_2 BLS_Period_invtransit_3_2 \
BLS_deltaChi2_invtransit_3_2 BLS_MeanMag_3_2 Mean_Mag_3 RMS_3 Expected_RMS_3 Npoints_3
----  ---------- ----- -------------- --------- ------------- --------- -------------- ---------- ---------- \
----------- ------------- ------------- ---------- ---------- ----------------- -------------------- \
----------------------- ----------------- --------------------------- --------------------------- \
---------------- ------------------ ------------------------- -------------- ---------- ---------- \
----------- ------------- ------------- ---------- ---------- ----------------- -------------------- \
----------------------- ----------------- --------------------------- --------------------------- \
---------------- ------------------ ------------------------- ------------------------- \
---------------------------- --------------- ---------- ----- -------------- ---------
EXAMPLES/3.transit 10.16727 0.00542 0.00104 3417 10.16714 0.00471 2.12298216 72.07052 0.00243 4.84830 \
0.01162 0.03000 0.98500 0.01000 -25264.50212 0.42585 146 4 93 87 0.00140 0.00407 14.99298 1.35250338 \
59.87855 0.00204 3.55079 0.00999 0.04500 0.54000 0.58000 -17824.72077 0.53476 132 4 130 66 0.00159 \
0.00437 11.36141 2.12704831 -2617.99005 10.16739 10.16664 0.00407 0.00104 3417

Run TFA and BLS on the light curves listed in EXAMPLES/lc_list_tfa (the only light curve in that list is EXAMPLES/3.transit). For TFA we use the trendlist given in EXAMPLES/trendlist_tfa (the second and third columns are the x and y positions for each star), we have TFA only include trend stars that are more than 25 pixels from the EXAMPLES/3.transit (the coordinates for this star are given in EXAMPLES/lc_list_tfa). We pass the corrected light curve to the subsequent commands and do not output the TFA coefficients or the TFA model. We then run BLS on the light curve cleaned by TFA. We search over a fixed q range of 0.01 to 0.1 and a fixed period range of 0.5 to 5.0 days. We search 2000 frequency points with 200 bins per period. We add +7 to the UT time to get the local time (used to determine the fraction of delta-chi2 that comes from one night), and report the top two peaks in the BLS spectrum. We do not output the periodogram, but do output the bls model light curve in the third column of EXAMPLES/OUTDIR1/3.transit.bls.model. We subtract the best fit box-car transit before passing the light curve to the next command. We also compute the rms before running TFA and after running BLS. The best period is in the column BLS_Period_1_2 (the _2 just denotes that this is the third command) and is 2.12298216 days in this case (with a delta-chi2 of -25264.50212 given in the BLS_deltaChi2_1_2 column), the second best period is 1.35250338 days. The SN for the first peak (72.07052) is the signal-to-noise measured in the BLS spectrum, the more commonly used signal to pink-noise ratio for the transit is 14.99298. Note that there are 146 points in transit, 4 observed transits and the estimated red and white noises of the light curve (after subtracting the model transit) are 0.00140 and 0.00407 mag respectively. The depth of the first peak is 0.01162 mag and the q is 0.03. Note that TFA reduces the RMS from 0.00542 mag (column RMS_0) to 0.00471 mag (column TFA_RMS_1), and subtracting the transit reduces it further to 0.00407 mag (column RMS_3).

$ ./vartools -i EXAMPLES/M3.V006.lc -Killharm fix 1 0.514333 10 0 1 \
   EXAMPLES/OUTDIR1/ fitonly outRphi -header
#Name Killharm_Mean_Mag_0 Killharm_Period_1_0 Killharm_Per1_Fundamental_Amp_2_0\
Killharm_Per1_Fundamental_Phi_2_0 Killharm_Per1_Harm_R_2_1_0 Killharm_Per1_Harm_Phi_2_1_0 \
Killharm_Per1_Harm_R_3_1_0 Killharm_Per1_Harm_Phi_3_1_0 Killharm_Per1_Harm_R_4_1_0 \
Killharm_Per1_Harm_Phi_4_1_0 Killharm_Per1_Harm_R_5_1_0 Killharm_Per1_Harm_Phi_5_1_0 \
Killharm_Per1_Harm_R_6_1_0 Killharm_Per1_Harm_Phi_6_1_0 Killharm_Per1_Harm_R_7_1_0 \
Killharm_Per1_Harm_Phi_7_1_0 Killharm_Per1_Harm_R_8_1_0 Killharm_Per1_Harm_Phi_8_1_0 \
Killharm_Per1_Harm_R_9_1_0 Killharm_Per1_Harm_Phi_9_1_0 Killharm_Per1_Harm_R_10_1_0 \
Killharm_Per1_Harm_Phi_10_1_0 Killharm_Per1_Harm_R_11_1_0 Killharm_Per1_Harm_Phi_11_1_0 \
Killharm_Per1_Amplitude_0
EXAMPLES/M3.V006.lc  15.77123     0.51433300   0.38041  -0.07662   0.47077   0.60826   0.35917   \
0.26249   0.23631  -0.06843   0.16353   0.60682   0.10621   0.28738   0.06203   0.95751   0.03602   \
0.58867   0.02900   0.22322   0.01750   0.94258   0.00768   0.66560   1.11128

$ echo EXAMPLES/4 | gawk '{amp = 0.25; for(i=1; i <= 10; i += 1) {print $1, amp; amp = amp*0.5;}}' \
     | ./vartools -l - -Injectharm fix 0.514333 10 \
        ampspec phaserand \
        ampfix 0.47077 amprel phasefix 0.60826 phaserel \
        ampfix 0.35916 amprel phasefix 0.26249 phaserel \
        ampfix 0.23631 amprel phasefix -0.06843 phaserel \
        ampfix 0.16353 amprel phasefix 0.60682 phaserel \
        ampfix 0.10621 amprel phasefix 0.28738 phaserel \
        ampfix 0.06203 amprel phasefix 0.95751 phaserel \
        ampfix 0.03602 amprel phasefix 0.58867 phaserel \
        ampfix 0.02900 amprel phasefix 0.22322 phaserel \
        ampfix 0.01750 amprel phasefix 0.94258 phaserel \
        ampfix 0.00768 amprel phasefix 0.66560 phaserel \
        0 0 \
        -LS 0.1 10.0 0.01 2 0 \
        -aov_harm 2 0.1 10.0 0.1 0.01 2 0 -header -numbercolumns

#1_Name 2_Injectharm_Period_0 3_Injectharm_Fundamental_Amp_0 4_Injectharm_Fundamental_Phase_0 \
5_Injectharm_Harm_2_Amp_0 6_Injectharm_Harm_2_Phase_0 7_Injectharm_Harm_3_Amp_0 8_Injectharm_Harm_3_Phase_0 \
9_Injectharm_Harm_4_Amp_0 10_Injectharm_Harm_4_Phase_0 11_Injectharm_Harm_5_Amp_0 12_Injectharm_Harm_5_Phase_0 \
13_Injectharm_Harm_6_Amp_0 14_Injectharm_Harm_6_Phase_0 15_Injectharm_Harm_7_Amp_0 16_Injectharm_Harm_7_Phase_0 \
17_Injectharm_Harm_8_Amp_0 18_Injectharm_Harm_8_Phase_0 19_Injectharm_Harm_9_Amp_0 20_Injectharm_Harm_9_Phase_0 \
21_Injectharm_Harm_10_Amp_0 22_Injectharm_Harm_10_Phase_0 23_Injectharm_Harm_11_Amp_0 24_Injectharm_Harm_11_Phase_0 \
25_LS_Period_1_1 26_Log10_LS_Prob_1_1 27_LS_SNR_1_1 28_LS_Period_2_1 \
29_Log10_LS_Prob_2_1 30_LS_SNR_2_1 31_Period_1_2 32_AOV_HARM_1_2 \
33_AOV_HARM_SNR_1_2 34_AOV_HARM_NEG_LOG_FAP_1_2 35_Period_2_2 36_AOV_HARM_2_2 \
37_AOV_HARM_SNR_2_2 38_AOV_HARM_NEG_LOG_FAP_2_2 39_Mean_lnAOV_2 40_RMS_lnAOV_2 \
EXAMPLES/4 0.51433300 0.25000 0.84019 0.11769 2.28864 0.08979 2.78305\
0.05908 3.29232 0.04088 4.80776 0.02655 5.32851 0.01551 6.83882\
0.00901 7.31017 0.00725 7.78491 0.00438 9.34446 0.00192 9.90766\
0.51493297 -492.06805 25.15612 1.06567664 -378.44559 19.24304 0.51432840 4322.49\
149.531 2969.03 0.51500100 3826.26 132.25 2805.09 28.5848 28.7159\
EXAMPLES/4 0.51433300 0.12500 0.39438 0.05885 1.39703 0.04489 1.44564\
0.02954 1.50910 0.02044 2.57873 0.01328 2.65368 0.00775 3.71819\
0.00450 3.74373 0.00363 3.77267 0.00219 4.88641 0.00096 5.00381\
0.51408199 -487.91541 25.18381 0.33926383 -379.53859 19.49478 0.51432840 4609.61\
137.471 3056.79 0.51414979 4585.38 136.744 3049.56 30.8468 33.3071\
...
...
EXAMPLES/4 0.51433300 0.00098 0.27777 0.00046 1.16381 0.00035 1.09581\
0.00023 1.04267 0.00016 1.99569 0.00010 1.95403 0.00006 2.90193\
0.00004 2.81087 0.00003 2.72319 0.00002 3.72033 0.00001 3.72112\
0.93527063 -69.42762 7.34581 0.33952304 -67.04912 7.08529 1.99136298 166.558\
9.40744 291.44 1.98761005 157.378 8.84087 276.201 14.1343 16.2025\
EXAMPLES/4 0.51433300 0.00049 0.55397 0.00023 1.71620 0.00018 1.92440\
0.00012 2.14745 0.00008 3.37667 0.00005 3.61120 0.00003 4.83530\
0.00002 5.02043 0.00001 5.20895 0.00001 6.48228 0.00000 6.75927\
0.93667874 -57.95590 7.76038 1.04106764 -55.69804 7.45058 0.99348763 199.742\
13.8977 345.37 0.99504450 188.754 13.0856 327.71 11.7027 13.5302\

In this example we conduct a simple test of injecting an RR Lyrae signal into a light curve with a fixed period and a range of amplitudes and attempt to recover the signal. We first fit a Fourier series to the light curve of an actual fundamental mode RR Lyrae star (M3 Variable V006). We do this using the -Killharm command fitting for the fundamental mode plus 10 harmonics. We fix the period to 0.514333 days which is the known period for this star. We use the keyword outRphi to output the fourier components as R_k1 = (A_k / A_1) and phi_k1 = (phi_k - k*phi_1) rather than the default sin and cos coefficients. The Fourier series in this format can then be used as a model signal to inject into a test light curve. This is done with the second call to vartools. The first call to echo/gawk generates an input list with the name of the input light curve (EXAMPLES/4) in the first column and the amplitude of the signal in the second column. This is piped into vartools (giving "-" as the argument to "-l" tells the program to read the list from stdin). We run three commands on each light curve: -injectharm, -LS and -aov_harm. The first command injects the model Fourier series into the light curve. We fix the period for this signal to 0.514333 days (alternatively we could have used "rand" 0.4 0.7 to inject random periods between 0.4 and 0.7 days, for example). We use 10 harmonics (plus the fundamental mode) in the model. The next line tells the program how to generate the amplitude/phase for the fundamental and the 10 lines after that tell it how to generate the amplitudes/phases of the higher harmonics. For the fundamental mode we use ampspec to read in the amplitude from the input list and phaserand to choose a random phase. For the 10 harmonics we use ampfix and phasefix to fix the R_k1 and phi_k1 values for each mode to the values that we found from fitting the Fourier series to the real RR Lyrae. Here we specify amprel to tell the program that we are inputing R_k1 rather than A_k and phaserel to tell the program that we are inputing phi_k1 rather than phi_k. If we didn't use the relative amplitudes and phases we would have to manually adjust the amplitudes and phases of each mode to keep the signal shape constant while varying the total amplitude and phase. After injecting the signal into the light curve, we attempt to recover the period of the injected signal using the Lomb-Scargle algorithm (-LS command) and the multi-harmonic AoV algorithm (-aov_harm command). The first two and last two lines of the output are shown. Note that the actual numbers that you find may be different depending on your random number generator (to have a non-repeatable set of simulated phases you would use the "-randseed time" option). Columns 3 through 24 of the output give the amplitude (A_k) and phase (phi_k) of each harmonic in the injected signal. The most important of these numbers in our case is the third column which gives the amplitude of the fundamental mode, which is what we varied in the input list. The top period recovered with L-S (column 25) ranges from 0.514 to 0.515 for the first 8 simulations (down to an amplitude of 1.95 mmag for the fundamental mode, corresponding to a peak-to-peak amplitude of 5.7 mmag for the full signal), for the two lowest amplitude simulations (0.98 and 0.49 mmag, corresponding to peak-to-peak amplitudes of 2.86 and 1.43 mmag) the signal is not recovered with LS. The same is true when using multi-harmonic AoV (column 31), but in this case the recovered period is generally closer to the input period (e.g. for the first simulation L-S finds 0.51493297 days while AoV gives 0.51432840), this is not surprising since the signal is not purely sinusoidal. 2 examples of the injected signal.

#!/bin/sh
star="$1"

if [ -n "$star" ] ; then

    vartools -i $star \
	-savelc \
	-clip 5.0 1 \
	-savelc \
	-LS 0.1 100. 0.1 3 0 clip 5. 1 \
	-aov 0.1 100. 0.1 0.01 1 0 clip 5. 1 \
	-aov_harm 1 0.1 100. 0.1 0.01 1 0 clip 5. 1 \
	-restorelc 1 \
	-clip 10.0 1 \
	-BLS q 0.01 0.1 0.1 20. 10000 200 7 2 0 0 0 \
	-restorelc 2 \
	-changeerror \
	-autocorrelation 0. 30. 0.1 EXAMPLES/OUTDIR1/ \
	-nobuffer

else


    PEXEC=${HOME}/bin/pexec
    SELF=$0
    LCLIST=EXAMPLES/lc_list

    vartools \
	-savelc \
	-clip 5.0 1 \
	-savelc \
	-LS 0.1 100. 0.1 3 0 clip 5. 1 \
	-aov 0.1 100. 0.1 0.01 1 0 clip 5. 1 \
	-aov_harm 1 0.1 100. 0.1 0.01 1 0 clip 5. 1 \
	-restorelc 1 \
	-clip 10.0 1 \
	-BLS q 0.01 0.1 0.1 20. 10000 200 7 2 0 0 0 \
	-restorelc 2 \
	-changeerror \
	-autocorrelation 0. 30. 0.1 EXAMPLES/OUTDIR1/ \
	-headeronly \
	-nobuffer \
	-numbercolumns
    
    $PEXEC -f $LCLIST -e star -R -t -z 19 -n 4 -c eval $SELF \$star
fi

This example is for a script that runs the -LS, -aov, -aov_harm, -BLS and -autocorrelation commands on a list of light curves. The script uses the pexec program by András Pál for parallelization. This script is included in the package as EXAMPLES/example6.sh. The if statement splits the script into two branches. The first branch runs vartools on a light curve if a light curve name is passed to the script as an argument. If no arguments are passed to the script, the second, initialization branch runs. In this branch vartools is called with no input light curve and with the "-headeronly" option which outputs the header for the command. After this, pexec is called. The -f option tells pexec to read the list of arguments (in this case light curve names) from EXAMPLES/lc_list. Each line read from this file will be assigned to the environment variable "star", then pexec will call this script passing $star as an argument. The pexec call will run 4 processes in parallel, each process will be run at nice +19. For the vartools command we first save the initial state of the light curve, then apply iterative 5-sigma clipping to the light curve. We save the 5-sigma clipped light curve. We then run the -LS, -aov and -aov_harm period finding algorithms. We then restore the light curve to its state before the 5-sigma clipping and apply 10-sigma clipping and run BLS (BLS would be sensitive to eclipses. To search for eclipses we would want to use a less aggressive sigma-clipping). After BLS we restore the light curve to its state after the 5-sigma clipping was applied and replace the errors in the light curve with the light curve RMS. Finally we run the -autocorrelation command on the light curve which will output the autocorrelation function to the EXAMPLES/OUTDIR1 directory (see for example EXAMPLES/OUTDIR1/2.autocorr which is periodic and has a first peak at 1.23 days).

1. Format the help sections to make them more legible, create manpages
2. Use automake for portability
3. Add examples for additional commands (e.g. transit injection/recovery) and TFA_SR
4. Purge non-useful commands (e.g. is -ensemblerescalesig useful and properly debugged?)
5. Improve the non-linear fitting for routines such as -MandelAgolTransit (consider incorporating MCMC)
6. Add the "fixcolumn" <colname | colnum> option to other commands that take parameters from previous commands (e.g. -Killharm)
7. Add a peak finder to -autocorrelation
8. Extend the -o command to allow a list of output filenames or a specified output filename
9. Add an option to print out the expected format of an input list given a command sequence
10. Add a function to do gawk-like manipulations of a light curve, or other arbitrary function manipulations.
11. Add fourier filtering routines.
12. Add a command to make plots of collections of light curves.
13. Add a routine to pick TFA template light curves, or generate simulated trend sets.
14. Light curve simulator