The dimension of the gamma ray arrival times array `TB' is set to be 10000 in the program, `TB(10000)'. This dimension can be easily changed by the user to accommodate their data.
Aside from the set of gamma ray arrival times, the user should also specify the range of the period to be searched by the program. Furthermore, the user should also specify the first and the second derivatives of the period, an epoch to be used as the reference point for the calculation of the phase value of each individual gamma ray arrival time as a fraction of one period, plus a housekeeping character string for the user to use. The assignment of an epoch value is optional. if a null value, i.e., a value less than 1 is assigned, the program will take the average value between the first and last data points to be the epoch. The input block also involves a multiplicative factor to be used to change the step size of the search effort. The use of this factor is explained in section 2.
While it may be more convenient for the user to express a periodicity structure in terms of its period, the mathematical analysis can be carried out more easily in terms of the frequency. After reading in the input data, the program calculates the corresponding frequency parameters from the period parameters for the subsequent analysis. Then the program determines the step size in frequency for the search effort in the way that if an error is made in deciding on the frequency by the amount of the step size, the accumulated error in the phase value from the earliest data point to the latest data point of the gamma ray arrival times is no bigger than one bin width in the phase plot (the light curve). This step size can be further modified with the multiplicative factor mentioned at the end of last paragraph. The program then determines the number of steps to be scanned over simply as the ratio of the search range to the step size. Afterwards, the main folding loop starts to cycle over the frequency steps.
The major part in the main folding loop is the determination of the phase values of the data points and the construction of the light curve. The program also calculates the Fourier power in the main folding loop. Because the overall observation time of a typical EGRET sighting is usually many times longer than the possible periods of interest, the calculation of the individual phase values must be done in double precisions to avoid cumulative errors from the rounding inaccuracies. Thus all the constants involved in the phase value calculations are declared to be REAL*8, including all the input data. On the other hand, once the phase value of a data point is determined, the calculation of its contribution to the Fourier power does not need such high precisions. Two arrays for the values of the sine and the cosine functions, 1000 elements each, in the angular range from 0 to 2p, are determined before the main folding loop starts. These arrays are then used in the place of the intrinsic functions themselves when sine or cosine values are needed in the Fourier power calculation. This standard practice to save computer CPU time when sine or cosine functions are called frequently is adopted in this program.
The output of this program involves no devices other than the line printer. Some preliminary information regarding this search effort, including the specified search range, the first and the second derivatives of the period, and the epoch, as well as some calculated values like the step size, number of steps, etc., are printed first. The housekeeping character string is also printed. Some simple information regarding the SELECT data file, including the number of records and the number of event records that have been read, and the header of the SELECT file, are written next. Then comes the main output block for all the vital information coming out of the search effort, one line for each frequency step. The information includes the numerical values of the light curve and the Fourier power, some indication of the significance of a possible periodicity exhibited in the light curve and a printer plot for this indication of significance, and a similar indication of significance and a printer plot for it for the Fourier power. The items in the main output block are arranged in a way that the first 72 columns contain the most interesting part for ready viewing on a remote terminal. Hard copies can be made from the output file at the user's discretion.
There is no way to foretell how much CPU time is needed for a specific search job. The user is advised to make some trial runs to see if the search effort is compatible with the computer time at the user's disposal.
The gamma ray arrival times come in through the SELECT data file. The connection to the SELECT data file is explained below for card no. 8.
Card no. 1 contains the epoch as the reference point for phase value calculations, in Julian Days, in format D17.10. The assignment of this epoch value is optional. If a null value is assigned (a value less than 1.0D0), the program will take the average value between the first and the last data points to be the epoch used in the program.
Card no. 2 contains the starting value (PERST) of the period search range, given in seconds, in format D17.10.
Card no. 3 contains the ending value (PEREND) of the period search range, given in seconds, in format D17.10.
Card no. 4 contains the first derivative of the period (PDOT), given in sec/sec, in format D17.10.
Card no. 5 contains the second derivative of the period (PDDT), given in sec/sec/sec, in format D17.10.
Card no. 6 contains the multiplicative factor (SDFRQ) used to modify the frequency step size as determined by the program, dimensionless, in format D17.10.
Now a word for this factor. The frequency step size is determined by the program in a unique way best suited for the period search effort. However, if the user wants to speed up the search or to refine the search, the step size as determined by the program can be modified with this multiplicative factor. This change in the step size will then be noted in the output file. However, if no such modification is desired, this factor should be set to unity, i.e., in the form +1.0000000000D+00.
Card no. 7 contains a string of 72 characters. This character string is provided for the user to note down any housekeeping items.
Card no. 8 contains the filename, in the form of a character string, of the SELECT data file from which the gamma ray arrival times, already corrected to the barycenter, are taken. This is the way the SELECT data file is connected to the program execution. This character string can be as long as necessary for the SELECT data filename.
Card no. 9 contains a single character to specify whether the SELECT data file is in the IBM format as when it was created by the IBM machine, or it has already been converted into the Sun format. A character `e' or `E' is for the IBM format (EBCDIC). A character `a' or `A' is for the Sun format (ASCII).
Afterwards, the main content of the output file for the result of this period search, one line for each search step, is printed out under a header of 9 items:
Item 1 is the period of this search step.
Item 2 is the chi square of the light curve, computed relative to the mean value of the count per bin.
Item 3 is the number of times the light curve crosses the mean value.
Item 4 is the probability that the null hypothesis that the light curve is a random distribution is true. This probability also takes into account the value of item 3. (When a peak exits in the light curve, i.e., when the data have a periodicity structure at this frequency, the number of times that the light curve crosses the mean value should be vastly reduced from a corresponding random distribution. The probability given in Item 4 takes this fact into consideration.)
Item 5 is a printer plot display for the probability of Item 4. What is plotted is the value of -log10(item 4). The height of this plot indicates the significance of a possible periodicity at this frequency.
Item 6 is the numerical value of the Fourier power.
Item 7 is an indication of the significance of a possible periodicity at this frequency. (See References.)
Item 8 is a printer plot display of the Fourier power. (See References.)
Item 9 is the numerical values of the light curve, divided into 20 bins.
While some meaning can be assigned to the significance of a possible periodicity at a certain period or frequency, the judgment of whether there is a real periodicity structure in the data or not must be left to the user. This programmer will not attempt to inject his own opinion into this question. Instead references are given below for the users to consult for a final decision on the result of the period search. The unpublished COS-B reports cited here are provided as part of the SEARCH program documentation set.
(1) The Source code consists of two files. One file is the main program plus a short subroutine to convert the EBCDIC character set into the ASCII character set. This file is called `search.f', and is written in Fortran77. The other file is a subroutine to convert the IBM floating point number representation into the Sun standard. This file is called `fpconv.c', and is written in C. This source code is stored at GSFC as part of the EGRET software collection.
(2) The input file consists of the items described in Section 2 above. It is the user's responsibility to supply this file. In fact, the content of this input file defines a particular periodicity search job.
(3) The output of the SEARCH program will be written on an output file. The user should assign a filename to receive the outcome of the program output.
(4) The SELECT data file provides the gamma ray arrival times to be processed by the SEARCH program. What SELECT data to use obviously depends on the user's own choice.
(5) The MAKEFILE will compile the source code, do the linking, and also define the command `search' to carry out the program execution. The MAKEFILE should look like this:
search: fpconv.o search.f
f77 -o search search.f fpconv.o
To execute the program, it is recommended to use the following command at the command line, after running the MAKEFILE:
search < [input filename] > [output filename]