Notes on export version of Program for Variable Star Period Search using Template Fitting.
This program fits the observed data for a varible star with 10 templates over range of potential periods, and locates most likely periods from the chi-squared minima. It produces data files and SuperMongo macros which help you visually select likely periods for further investigation.
I find this program useful when (a) you have some external evidence of what type of variables are in your data set (e.g., you are looking at the mag/color of the RR Lyrae in a globular cluster); and (b) your data are undersampled, so less biased techniques like Phase Dispersion Minimization, String Length, or Lomb-Scargle Method fail.
For more information, see the relevant sections of Layden, Ritter,Welch, & Webb (1999, AJ, 117, 1313) and Layden & Sarajedini (2000, AJ, 119, 1760).
NB: compile w/ -e option:
%f77 stdlc_per10_export.f -e -o stdlc_per10_export.e
% stdlc_per10_export.e
* READS template file (stdlc.all10) and reports
* Requests which observing runs to use. This is useful because it is sometimes easiest to isolate an approximate period using one or two closely-spaced runs of data, and then do higher resolution period searches around the likely periods using more/all of the runs. Reads info on runs from stdlc_per10.runs (which you modify to suit your observational data set). Here, you indicate which runs you want to include in the analysis (1) and which you do not (0).
Number of observing runs = 5 Run(i), Time(begin), Time(end) = 1 520.2 522.2 Use data from this run? (1=yes, 0=no): 1 Run(i), Time(begin), Time(end) = 2 529.2 531.2 Use data from this run? (1=yes, 0=no): 1 Run(i), Time(begin), Time(end) = 3 638.2 642.2 Use data from this run? (1=yes, 0=no): 1 Run(i), Time(begin), Time(end) = 4 896.2 898.2 Use data from this run? (1=yes, 0=no): 1 Run(i), Time(begin), Time(end) = 5 1278.2 1279.2 Use data from this run? (1=yes, 0=no): 0
* Requests name of data containing observed data:
Enter name of file w/ HJD Vmag Imag (vfe.B01): vfe.A01b
* Requests "weighting factor" for data points taken under conditions of bad seeing. For example, if you specify 2.0, it will multiply by 2.0 the error estimates of all the Vmags which you have designated as bad seeing (see stdlc_export.f), so they count only 1/2^2 as much as the data taken during conditions of good seeing. This is useful if you are working in a crowded field where stellar profiles become blended in bad seeing.
Enter error ratio for bad seeing points (eg. 2.0): 1.5
* (New addition 2005): Requests "maximum error to include". The program will exclude any magnitudes in the first column of the input file that have photometric errors larger than this value (there is also a "backstop" value of 1.0 to eliminate any missing ("out of bounds") magnitude/error values in the third/fourth columns of the input file).
Enter maximum error value to accept (eg. 0.1): 0.05
* Reports how many data points were read:
READ 12 mags from variable star data file
* Requests name for output file:
Enter name for output file: sample.A01b
* Requests how much detail you want reported during fitting process.
Output verbosity in reporting indiv solns (hi/lo/no = 2/1/0): 0
* Enter the period where you want to begin fitting (eg 0.2 days), where you will end fitting (eg, 1.0 days), and the fractional increment (ie, compute the next period, P1, from the current period, P0, according to P1={1.0+INCR}*P0). Obviously, a smaller value for the increment will take longer to run, but provide higher resolution results. You must determine what the best increment is for your data set.
Enter Starting and Ending Periods, and Fractional Increment (eg 0.2 1.0 0.001): 0.2 1.0 0.005
* THINK THINK THINK. Starting at the beginning period (eg. 0.2 days), it fits each of the 10 templates to the observed data, and computes a value of chi-squared for each fit. It then determines the next period to try and does the fitting there. Ad nauseum. It really is a big hammer, pounding away at a small problem, but heck, that's what we have computers for! It spits a progress report at you every 300th period:
i per(i) = 300 0.2699351 i per(i) = 600 0.3643245 i per(i) = 900 0.4917198 i per(i) = 1200 0.6636626 i per(i) = 1500 0.8957289 Ending simulation, see sample.A01b
* Reports format of outupt file, in this example, sample.A01b, and the number of periods it examined. Also tells the largest chi-squared value it encountered and suggests an upper limit for the plot window (old stuff, but useful if you will be plotting manually).
Output = Period xsq(1) xsq(2) xsq(3) xsq(4) xsq(5) xsq(6) xsq(7) xsq(8) xsq(9) xsq(10) N Per Pts = 1610 Maximum chi^2 and suggested plot limit = 0.3109 0.1891
* Reminds you what commands to type if you are using SuperMongo to plot chi-squared as a function of period for each template. The program has created 2 files, in this example (1) sample.A01b (see above) and (2) sm.sample.A01b, a SM macro that plots the result.
To plot with SM: macro read sm.sample.A01b plt
* Upon viewing this plot, you can look for minima (ie, periods that produced good fits). Occasional upward or downward spikes result when the fitting failed to converge on a realistic solution, so generally should be ignored. Sometimes, there is a "packet" of minima, presumably where data in one or more closely-spaced observing runs "beat" with data in one or more distantly-spaced observing runs. Here is a plot showing the result for this example (encapsulated postscript).
In our example, the best period appears to be obtained with Template#2, at Period ~ 0.5 days.
* I have several short programs which you may run to help you pick out the minima (ie, most likely period). See below.
* Once you have determined the most likely period, use stdlc_export.f to plot the final light curve. It will plot the data (both the"primary" magnitudes in the third column of the input data file and the "secondary" magnitudes from other filters that may appear in the fifth and seventh columns of the input data file) along with the best-fittingtemplate.
* In general, I have found that the fit with the smallest chi-squared value among all of period-template space does in fact produce the cleanest appearing light curve. It also tends to find the true period of a star whose light curve I have deliberately undersampled (eg, used only a random sampling of 1/10 of the available data points). I therefore believe the program to be good period-finding tool. However, I have not done more elaborate statistical simulations. You may want to do that on your data sets. Also, you should do some careful tests on your own or simulated data to ensure you fully understand what the program is doing. You use this program at your own risk.
* Acknowledgement: If you use this program in published work, I would appreciate your citing it as "Layden, A.C., Ritter, L.A., Welch, D.L., & Webb, T.M.A. 1999, AJ, 117, 1313"
Once again, please email me if you have questions, comments, or problems.
Best of luck,
Andy Layden
ACL -- updated 2005 Apr 14
After viewing the plot, you decide you want to find the exact period corresponding to particular minimum. Specify the lower and upper bounds of the period range you wish to search. Program returns the minimum chi-squared value (chimin) in each of the 10 templates (col) and the period at which it occurs, per(chimin).
* Enter name of file containing chi-squared(period) data:
Enter name of STDLC_PER output: sample.A01b Enter range of periods to consider (Pmin Pmax): 0.4 0.6 Col chimin per(chimin) = 1 0.0015770 0.5026525 Col chimin per(chimin) = 2 0.0062000 0.5011474 Col chimin per(chimin) = 3 0.0043560 0.5021503 Col chimin per(chimin) = 4 0.0130740 0.5001466 Col chimin per(chimin) = 5 0.0144690 0.5001466 Col chimin per(chimin) = 6 0.0219710 0.5001466 Col chimin per(chimin) = 7 0.0452260 0.5071947 Col chimin per(chimin) = 8 0.0541050 0.4936897 Col chimin per(chimin) = 9 0.1473260 0.5107559 Col chimin per(chimin) = 10 0.1296460 0.4431711
In this example, the template with the lowest chi-squre value is #1 (col=1), where chimin=0.0015770 at a period of 0.5026525 days. It is a good candidate period for the variable star. You may want to repeat this for several different minima, and compare the appearances of their light curves at these fixed periods using stdlc_export. I always visually inspect each light curve this way, to ensure I am not misled by a false-positive chi-square value.
For a given template, which points fall below some maximum value of chi-squared (which you specify)? This can be useful in determining which is the smallest chimin in a "packet", as discussed above.
Enter name of STDLC_PER output: sample.A01b Which x^2 column should I consider (1,2,...,ncol): 3 Find minima below this value of x^2 (eg 0.05): 0.01 chimin per(chimin) = 0.0098390 0.4941834 chimin per(chimin) = 0.0086700 0.4951723 chimin per(chimin) = 0.0076470 0.4961632 chimin per(chimin) = 0.0067560 0.4971561
...etc...
In this example, we find all the periods for which chimin is less than 0.01 using template #3.