Pulsar Tools Anatomy
This article describes data processing procedures by the pulsar analysis tools distributed as a part of a Science Tools package. The versions of the tools described here is the ones included in Science Tools version v9r8p2: gtbary version v6r2p4, gtephem version v8r2, gtpulsardb version v8r2, gtpphase version v8r3, gtophase version v8r3, gtpsearch version v10r5, gtpspec version v10r5, gtptest version v10r5. A different version of those tools may or may not behave the same way as described in this article. All of the above tools report their version numbers on the first line of the tools' textual output. If you find the tool behave differently than described in this article, please check the version number of the tool in its textual output against the ones listed above.
Chapter 1. Data processing procedures of pulsar analysis tools
In this chapter, data processing procedures are described step by step in great detail for gtephem, gtpphase, gtophase, gtpsearch, gtpspec, and gtptest. Those tools are designed to play central roles in users' pulsar analysis. In this chapter, data processing procedures by those tools are described in detail step by step.
gtephem
 The tool determines the name of a leap second file to use in the later stage of the process when necessary as described in Step 2 of Chapter 3.
 The tool loads pulsar ephemerides as described in Step 4 of Chapter 3.
 The tool computes a reference time of computation, or a moment in time that reftime parameter points to as in a time format specified by timeformat parameter in a time system specified by timesys parameter.
 The tool chooses the best spin ephemeris for the reference time as described in Appendix C. If strict parameter is set to true, ephemeris validity is required. Otherwise, ephemeris validity is not required.
 The tool displays the chosen spin ephemeris.
 The tool chooses the best orbital ephemeris for the reference time as described in Appendix C.
 The tool displays the chosen orbital ephemeris.
 With the chosen spin ephemeris, the tool computes and displays the following quantities at the reference time, as described in Appendix A: the pulsar's Right Ascension and Declination, a pulse phase, a pulse frequency, the first timederivative of pulse frequency, and the second timederivative of pulse frequency.
 The tool searches for an ephemeris remark that is effective between the reference time and the reference epoch of the chosen ephemeris. If found, the tool produces a warning message listing ephemeris remarks that are found.
gtpphase
 The tool opens input event file(s) as described in Step 1 of Chapter 3.
 The tool determines the name of a leap second file to use in the later stage of the process when necessary as described in Step 2 of Chapter 3.
 The tool interprets user's request on arrival time corrections as described in Step 3 of Chapter 3.
 The tool loads pulsar ephemerides as described in Step 4 of Chapter 3.
 The tool initializes arrival time corrections as described in Step 5 of Chapter 3.
 The tool summarizes ephemeris status as described in Step 6 of Chapter 3.
 The tool reads pphasefield parameter from the parameter interface of the tool, and use it for the name of a FITS column which the tool writes pulse phase values.
 For all given event files, the tool creates a new FITS column with the name if it doesn't exist.
 The tool reads pphaseoffset parameter from the parameter interface of the tool, and keeps the value during the process.
 For each event in the given event file(s), the tool takes the following steps to assign a pulse phase to the event.
 The tool reads a value from a FITS column specified by timefield parameter, and computes a photon arrival time of the event.
 The tool applies arrival time corrections to the photon arrival time as described in Appendix D.
 The tool chooses the best spin ephemeris for the corrected arrival time as described in Appendix C requiring ephemeris validity.
 With the chosen spin ephemeris, the tool computes a pulse phase at the corrected time as described in Appendix A, adds pphaseoffset parameter value to it, then truncates the integer part of the value so that the result is nonnegative and less than one (1).
 The tool writes the truncated phase value to a FITS column specified by pphasefield parameter.
gtophase
 The tool opens input event file(s) as described in Step 1 of Chapter 3.
 The tool determines the name of a leap second file to use in the later stage of the process when necessary as described in Step 2 of Chapter 3.
 The tool interprets user's request on arrival time corrections as described in Step 3 of Chapter 3.
 The tool loads pulsar ephemerides as described in Step 4 of Chapter 3.
 The tool initializes arrival time corrections as described in Step 5 of Chapter 3.
 The tool summarizes ephemeris status as described in Step 6 of Chapter 3.
 The tool reads ophasefield parameter from the parameter interface of the tool, and use it for the name of a FITS column which the tool writes orbital phase values.
 For all given event files, the tool creates a new FITS column with the name if it doesn't exist.
 The tool reads ophaseoffset parameter from the parameter interface of the tool, and keeps the value during the process.
 For each event in the given event file(s), the tool takes the following steps to assign an orbital phase to the event.
 The tool reads a value from a FITS column specified by timefield parameter, and computes a photon arrival time of the event.
 The tool applies arrival time corrections to the photon arrival time as described in Appendix D.
 The tool chooses the best orbital ephemeris for the corrected arrival time as described in Appendix C.
 With the chosen orbital ephemeris, the tool computes an orbital phase at the corrected time as described in Appendix A, adds ophaseoffset parameter value to it, then truncates the integer part of the value so that the result is nonnegative and less than one (1).
 The tool writes the truncated phase value to a FITS column specified by pphasefield parameter.
gtpsearch
 The tool opens input event file(s) as described in Step 1 of Chapter 3.
 The tool determines the name of a leap second file to use in the later stage of the process when necessary as described in Step 2 of Chapter 3.
 The tool interprets user's request on arrival time corrections as described in Step 3 of Chapter 3.
 The tool loads pulsar ephemerides as described in Step 4 of Chapter 3.
 The tool initializes arrival time corrections as described in Step 5 of Chapter 3.
 The tool summarizes ephemeris status as described in Step 6 of Chapter 3.
 The tool chooses the best spin ephemeris for the time origin computed in Step 52 of Chapter 3.
 With the chosen spin ephemeris, the tool computes a pulse frequency at the time origin.
 The tool computes the Fourier resolution of the given event data in the following procedure.
 The tool computes the start time and the stop time of the given event data as described in Appendix B.
 The tool applies arrival time corrections to the start time and the stop time as described in Appendix D. If it is not successful, then the tool produces an error message and terminates the process.
 The tool computes a time difference between the corrected start time and the corrected stop time in the time system that is determined to be used for time series analysis in Step 52 of Chapter 3.
 The tool computes an inverse of the time difference computed in the previous step, and takes it as the Fourier resolution of the given event data.
 The tool determines trial frequencies in the following procedure.
 The tool reads numtrials parameter from the parameter interface of the tool, and sets it to the number of trial frequencies.
 The tool uses the pulse frequency at the time origin computed above to determine the center of trial frequency range. If the number of trial frequencies is odd, the computed frequency is set to the center of trial frequency range. If the number of trial frequencies is even, the computed frequency is set to the first frequency after the center of trial frequencies. For example, if the number of trial frequencies is 11, the computed frequency is set to the 6th trial frequency. If the number of trial frequencies is 12, the computed frequency is set to the 7th trial frequency.
 The tool reads scanstep parameter from the parameter interface of the tool, and multiplies it to the Fourier resolution computed in the previous step, and uses it as a frequency step size, or a difference in frequency between neighboring trial frequencies.
 The tool computes and sets all trial frequencies based on the central frequency and the frequency step size, both of which are computed above.
 The tool evaluates algorithm parameter, chooses a statistical test to perform on the given data, and initializes the chosen statistical test as shown in the following table.
algorithm value 
Statistical test 
Initialization 
CHI2 
χ^{2} test 
The tool reads numphase parameter from the parameter interface of the tool, and uses it for the number of phase bins of a folded light curve. 
RAYLEIGH 
Rayleigh test 
Same as Z_{n}^{2} test for n = 1. 
Z2N 
Z_{n}^{2} test 
The tool reads numharm parameter from the parameter interface of the tool, and uses it for the number of harmonics to sum up (i.e., the value of the suffix "n" of Z_{n}^{2} test). 
H 
H test 
The tool reads maxharm parameter from the parameter interface of the tool, and uses it for the maximum number of harmonics to evaluate in searching for the optimal number of harmonics to sum up. 
Those statistical tests are wellknown and their definitions and descriptions are found in the literature. For example, Leahy et al. (1983) described χ^{2} test in detail, and Vaughan et al. (1994) pointed out importance of noisesignal interactions in χ^{2} test. Also, Buccheri et al. (1983) explained Z_{n}^{2} test as a part of their data analysis, and De Jager et al. (1989) defined and discussed H test.
References
Buccheri et al. 1983, A&A 128, 245
De Jager et al. 1989, A&A 221, 180
Leahy et al. 1983, ApJ 266, 160
Vaughan et al. 1994, ApJ 435, 362
 For each event in the given event file(s), the tool takes the following steps to accumulate contribution(s) of the event to the test statistic for the chosen statistical test at all the trial frequencies.
 The tool reads a value from a FITS column specified by timefield parameter, and computes a photon arrival time of the event.
 The tool applies arrival time corrections to the photon arrival time as described in Appendix D.
 The tool computes an elapsed time at the corrected arrival time since the time origin computed in Step 52 of Chapter 3. The computation is performed in a time system that is determined to be used for time series analysis in Step 52 of Chapter 3.
 For each of all the trial frequencies, the tool computes a contribution of this event to the test statistic of the chosen statistical test. First, a pulse phase is computed for each trial frequency. Then, if χ^{2} test was chosen in the previous step, one event is added to an appropriate phase bin of a folded light curve for the trial frequency. If Z_{n}^{2} test, Rayleigh test, or H test was chosen, the tool uses the computed pulse phase to compute the sine and the cosine components of Z_{n}^{2} test for all harmonics that are needed to compute the chosen test statistic, and add them to the sum of those components over events that have already been processed.
 The tool computes values of the test statistic for the chosen statistical test for all trial frequencies.
 The tool identifies the trial frequency that gives the maximum statistic, and evaluates a chance probability that a value of the test statistic of the chosen statistical test exceeds the maximum value found under the given test condition, under the assumption that the no pulsations exist in the given time series data. The procedure to compute it can be further broken down as follows.
 The tool computes a singletrial chance probability, p, which is a probability that a value of the test statistic of the chosen statistical test exceeds the maximum statistic identified in the previous step after one trial. If no pulsations exist in the given time series data, the test statistic follows a probability distribution shown in the table below, where numphase and numharm are values of the named parameters taken from the parameter interface of the tool.
algorithm value 
The test statistic follows... 
CHI2 
χ^{2} distribution with (numphase  1) degrees of freedom 
RAYLEIGH 
χ^{2} distribution with 2 degrees of freedom 
Z2N 
χ^{2} distribution with numharm×2 degrees of freedom 
H 
special probability distribution for H test, defined by De Jager et al. (1989) 
For a technical reason, the tool computes a range of probability that the true value of the singletrial chance probability lies. The algorithm to compute a probability range for a given statistic value is described in a separate document titled "Numerical integration of the chisquared probability density function" (PDF file, 148 kB) for your reference.
References
De Jager et al. 1989, A&A 221, 180
 The tool computes the number of independent trials, N, which is the number of trial frequencies that are considered statistically independent of each other. Letting f_{min} be the smallest trial frequency, f_{max} the largest, and f_{Fourier} the Fourier resolution of the given event data, the tool first computes the maximum number of independent trial frequencies that can fit in the range of trial frequency, N_{max} , by:
N_{max} = ceil((f_{max}  f_{min}) / f_{Fourier})
where ceil(x) is the ceiling function of x whose value is the smallest integer that is greater than or equal to x. Then, the tool compares it with the actual number of trial frequencies, N_{trial} (which is specified by numtrials parameter in the parameter interface of the tool), and determines N to be the smaller of them, or:
N = min(N_{max} , N_{trial})
 The tool computes a multitrial chance probability, Q, which is defined as:
Q = 1  (1  p)^{N}
whose computation procedure is described in a separate document titled "Numerical computation of false alarm probability for multiple trials" (PDF file, 148 kB) for your reference. Note that Q is computed for each of the lower bound of p and the upper bound, and the tool takes them as the lower and the upper bounds of Q, respectively.
 The tool evaluates outfile parameter, and creates an output FITS file unless the value of the parameter is NONE. The value of outfile parameter is used for the name of the output file, and the output file contains the test result, including values of the test statistic at all trial frequencies, a multitrial chance probability computed in the previous step, and other information related to the statistical test.
 The tool reports the same test result as in an output FITS file for human reading.
 The tool evaluates plot parameter and title parameter to create a plot showing the result.
 If the value of plot parameter is a logical true, the tool creates a graphical output that plots values of the test statistic against trial frequencies. Otherwise, the tool does nothing in this step.
 If the value of title parameter is "DEFAULT", the tool uses a predefined title for a title of the plot. Otherwise, the value of title parameter is used for the plot title.
gtpspec
 The tool opens input event file(s) as described in Step 1 of Chapter 3.
 The tool determines the name of a leap second file to use in the later stage of the process when necessary as described in Step 2 of Chapter 3.
 The tool interprets user's request on arrival time corrections as described in Step 3 of Chapter 3.
 The tool loads pulsar ephemerides as described in Step 4 of Chapter 3.
 The tool initializes arrival time corrections as described in Step 5 of Chapter 3.
 The tool summarizes ephemeris status as described in Step 6 of Chapter 3.
 The tool computes the start time and the stop time of the given event data as described in Appendix B.
 The tool applies arrival time corrections to the start time and the stop time as described in Appendix D. If it is not successful, then the tool produces an error message and terminates the process.
 The tool reads binwidth parameter from the parameter interface of the tool to define time bins for DFFT (Discrete Fast Fourier Transform) algorithm. A time bin is a time interval that the tool counts the number of photon arrival times that fall in. In this step, the tool defines time bins with the length specified by binwidth parameter, starting at the corrected start time.
 The tool reads numbins parameter from the parameter interface of the tool, and defines time segments for DFFT computations. A time segment is a group of successive time bins that DFFT is performed at once. In this step, the tool defines time segments such that the number of time bins in one time segment is equal to the value of numbins parameter, starting at the corrected start time.
 For each event in the given event file(s), the tool takes the following steps to determine a time segment and a time bin that the event belongs to.
 The tool reads a value from a FITS column specified by timefield parameter, and computes a photon arrival time of the event.
 The tool applies arrival time corrections to the photon arrival time as described in Appendix D.
 The tool computes an elapsed time at the corrected arrival time since the corrected start time. The computation is performed in a time system that is determined to be used for time series analysis in Step 52 of Chapter 3. Then:
 if the corrected photon arrival time is at or later than the corrected start time and at or earlier than the corrected stop time, the tool divides the computed elapsed time by the length of a time segment (= binwidth×numbins), and takes the ratio as the segment number and the reminder the time bin number.
 Otherwise, the event is discarded, i.e., the tool stops processing this event and moves on to the next event.
 For each time segment, the tool takes the following steps to compute a power spectral density using DFFT algorithm.
 The tool counts the number of events that belongs to each time bin in a time segment being processed to create a light curve, or an array whose element stores the number of events in a time bin.
 The tool computes the average number of events per time bin, and subtract the result from the element of the light curve computed in the previous step.
 The tool performs DFFT on the time segment.
 At every frequency of the computed Fourier transform of the light curve, the tool computes a Fourier power, or a squared sum of the real part and the imaginary part of the Fourier transform.
 The tool adds the Fourier power spectra for all the time segments computed in the previous step.
 The tool reads lowfcut parameter, and uses it for the lower boundary of frequency for a pulsation search. The tool identifies the lowest frequency bin in the computed Fourier transform whose frequency is not less than the lowfcut parameter value. In all the procedure described below, the tool ignores all frequency bins below the lowest frequency bin.
 The tool identifies the trial frequency that gives the maximum statistic, and evaluates a chance probability that a summed Fourier power exceeds the maximum value found under the given test condition. The procedure to compute it can be further broken down as follows.
 The tool computes a singletrial chance probability, p, which is a probability that a summed Fourier power exceeds the maximum Fourier power identified in the previous step after one trial. If no pulsations exist in the given time series data, a summed Fourier power follows the χ^{2} distribution with 2×M degrees of freedom, where M is the number of time segments that contains at least one event. For a technical reason, the tool computes a range of probability that the true value of the singletrial chance probability lies. The algorithm to compute a probability range for a given statistic value is described in a separate document titled "Numerical integration of the chisquared probability density function" (PDF file, 148 kB) for your reference.
 The tool counts the number of frequency bins that are searched for the maximum Fourier power, N, starting at the lowest frequency bin identified in the earlier step and ending at the last frequency bin which corresponds to the Nyquist frequency.
 The tool computes a multitrial chance probability, Q, which is defined as:
Q = 1  (1  p)^{N}
whose computation procedure is described in a separate document titled "Numerical computation of false alarm probability for multiple trials" (PDF file, 148 kB) for your reference. Note that Q is computed for each of the lower bound of p and the upper bound, and the tool takes them as the lower and the upper bounds of Q, respectively.
 The tool evaluates outfile parameter, and creates an output FITS file unless the value of the parameter is NONE. The value of outfile parameter is used for the name of the output file, and the output file contains the test result, including the summed Fourier powers at and above the lowest frequency determined in the earlier step, a multitrial chance probability computed in the previous step, and other information related to the pulsation search.
 The tool reports the same test result as in an output FITS file for human reading.
 The tool evaluates plot parameter and title parameter to create a plot showing the result.
 If the value of plot parameter is a logical true, the tool creates a graphical output that plots the summed Fourier powers against pulse frequencies. Otherwise, the tool does nothing in this step.
 If the value of title parameter is "DEFAULT", the tool uses a predefined title for a title of the plot. Otherwise, the value of title parameter is used for the plot title.
gtptest
 The tool opens input event file(s) as described in Step 1 of Chapter 3.
 The tool checks whether all the given event files contain a FITS column with the name specified by pphasefield parameter. If it is missing in one of the event files, the tool produces an error message and terminates the process.
 The tool prepares and initializes all the statistical tests that are currently implemented: χ^{2} test, Rayleigh test, Z_{n}^{2} test, and H test.
 The tool reads numphase parameter from the parameter interface of the tool, and the parameter value is used for the number of phase bins of a folded light curve in χ^{2} test.
 The tool reads numharm parameter from the parameter interface of the tool, and the parameter value is used for the number of harmonics to sum up in Z_{n}^{2} test (i.e., the value of the suffix "n" of Z_{n}^{2} test).
 The tool reads maxharm parameter parameter from the parameter interface of the tool, and the parameter value is used for the maximum number of harmonics to evaluate in searching for the optimal number of harmonics to sum up in H test.
 For each event in the given event file(s), the tool takes the following steps to accumulate contribution(s) of the event to the test statistics for all the statistical tests.
 The tool reads a FITS column specified by pphasefield parameter.
 The tool computes a contribution of this event to each test statistic of all the statistical test. For
χ^{2}, one event is added to an appropriate phase bin of a folded light curve. For Z_{n}^{2} test, Rayleigh test, and H test, the tool uses the pulse phase to compute the sine and the cosine components of
Z_{n}^{2} test for all harmonics that are needed to compute the test statistics, and add them to the sum of those components over events that have already been processed.
 The tool computes values of all the test statistics.
 The tool evaluates a chance probability that a value of each test statistic exceeds the value computed from the given event file(s), under the assumption that the no pulsations exist in the given time series data. Specifically speaking, null hypotheses for all the implemented statistical test are following.
 For χ^{2} test, the test statistic is assumed a sample from the χ^{2} distribution with
(numphase  1) degrees of freedom, where numphase is the value of numphase parameter taken from the parameter interface of the tool.
 For Rayleigh test, the test statistic is assumed a sample from the χ^{2} distribution with 2 degrees of freedom, where numharm is the value of numphase parameter taken from the parameter interface of the tool.
 For Z_{n}^{2} test, the test statistic is assumed a sample from the χ^{2} distribution with
(numharm × 2) degrees of freedom, where numharm is the value of numphase parameter taken from the parameter interface of the tool.
 For H test, the test statistic is assumed a sample from the special probability distribution for H test, defined by De Jager et al. (1989).
For a technical reason, the tool computes a range of probability that the true value of the singletrial chance probability lies. The algorithm to compute a probability range for a given statistic value is described in a separate document titled "Numerical integration of the chisquared probability density function" (PDF file, 148 kB) for your reference.
References
De Jager et al. 1989, A&A 221, 180
 The tool evaluates outfile parameter, and creates an output FITS file unless the value of the parameter is NONE. The value of outfile parameter is used for the name of the output file, and the output file contains the test result, including the folded light curve created for χ^{2} test, values of the test statistics, and other information related to the statistical tests.
 The tool reports the same test result as in an output FITS file for human reading.
 The tool evaluates plot parameter and title parameter to create a plot showing the result.
 If the value of plot parameter is a logical true, the tool creates a graphical output that plots the folded light curve created for χ^{2} test. Otherwise, the tool does nothing in this step.
 If the value of title parameter is "DEFAULT", the tool uses a predefined title for a title of the plot. Otherwise, the value of title parameter is used for the plot title.
Chapter 2. Other supporting tools for pulsar analysis
There are two more pulsarrelated tools in a Science Tools package: gtbary and gtpulsardb. Those tools use the same components for fundamental computations such as arrival time corrections, time system conversions, and pulsar ephemerides loading, as the ones used by the tools described in the previous chapter. But, due to their distinct functionality, their procedural steps are different. In this chapter, data processing procedures by the two supporting tools are described in detail step by step.
gtbary
 The tool evaluates tcorrect parameter to determine a type of arrival time corrections to perform. If tcorrect parameter is set to BARY, then the tool performs barycentric corrections. If tcorrect parameter is set to GEO, then the tool performs geocentric corrections.
 The tool opens an input event file specified by evfile parameter. It is not allowed to give more than one event file at a time.
 The tool copies the input event file to a temporary output file. The name of the temporary output file is the value of
outfile parameter with ".tmp" attached at the end. For example, if outfile parameter is set to
"gtbary_output.fits", then the temporary output file name is "gtbary_output.fits.tmp".
 The tool reads solareph parameter from the parameter interface of the tool, and determines a solar system ephemeris to use in barycentric corrections.
 The tool determines the name of a leap second file to use in the later stage of the process when necessary as described in Step 2 of Chapter 3.
 The tool reads ra and dec parameters from the parameter interface of the tool, and uses them as the pulsar's Right Ascension and Declination for which geocentric or
barycentric corrections are performed. Note that the values of ra and dec parameters are interpreted in the same celestial reference frame as the one used by the solar system ephemeris chosen in the earlier step. To be specific, FK5 is used by JPL DE200, and ICRS by JPL DE405.
 The tool updates values of the following header keywords of all extensions in the temporary output file.
Header keyword 
The value is set to... 
TIMESYS 
TDB if tcorrect = BARY, and TT if tcorrect = GEO. 
TIMEREF 
SOLARSYSTEM if tcorrect = BARY, and GEOCENTRIC if tcorrect = GEO. 
RA_NOM 
Right Ascension used for geocentric or barycentric corrections. 
DEC_NOM 
Declination used for geocentric or barycentric corrections. 
RADECSYS 
the name of celestial reference frame (FK5 or ICRS). 
PLEPHEM 
solar system ephemeris used for geocentric or barycentric corrections (JPLDE200 or JPLDE405). 
TIMEZERO 
0. (zero) 
CREATOR 
gtbary v6r2p4 
DATE 
file creation date in YYYYMMDDThh:mm:ss format. 
TIERRELA 
the same value as in the input event file, and 1.e9 if it is
missing in the input event file. 
 The tool reads scfile, sctable, and angtol parameters from the parameter interface of the tool, and stores them in a memory space for later use.
 The tool reads the following time expressions in the input event file, applies geocentric or barycentric corrections to them, and writes the result to the temporal output file.
 TSTART header keyword value
 TSTOP header keyword value
 DATEOBS header keyword value
 DATEEND header keyword value
 FITS column specified by timefield parameter of all extensions other than GTI extensions
 START column of GTI extensions
 STOP column of GTI extensions
Barycentric times are computed as described in Appendix D. Geocentric times are computed in the same way as barycentric times, except that a geocentric time, t_{geo} defined by the following formula, is computed in place of t_{bary} in Appendix D.
t_{geo} = t_{obs} + r_{sat} · n_{psr} / c
where r_{sat} is the vector pointing from the center of the Earth to the spacecraft. See Appendix D for definitions of the other quantities in the above formula.
 The tool renames the temporal output file to a file name specified by outfile parameter.
gtpulsardb
 The tool determines the name of a leap second file to use in the later stage of the process when necessary as described in Step 2 of Chapter 3.
 The tool evaluates psrdbfile parameter, and creates a list of file names from the parameter value. If the value of psrdbfile parameter starts with an atmark ("@"), the part that follows the atmark is interpreted as the name of a text file which contains a list of pulsar ephemerides database file names, and the contents of the text file is set to the file name list. Otherwise, the value of psrdbfile parameter is interpreted as a pulsar ephemerides database file name, and set to the first and only element of the file name list.
 The tool loads the entire contents of all the pulsar ephemerides database files in the file name list, in the order of their appearance in the file name list.
 The tool evaluates filter parameter, and filter the loaded pulsar ephemerides as follows. Then, the tool removes entries in OBSERVERS extension and ALTERNATIVE_NAMES extension that become unnecessary as a result of the filtering.
 If filter = NAME, the tool reads psrname parameter from the parameter interface of the tool, and filters the loaded pulsar ephemeris data by pulsar name. Spin ephemerides, orbital ephemerides, and ephemeris remarks for the pulsar specified by psrname parameter are subselected, and those for other pulsars are removed.
 If filter = TIME, the tool reads tstart and tstop parameters from the parameter interface of the tool, and filters the loaded pulsar ephemeris data by time. Spin ephemerides whose validity window overlaps with the time interval specified by tstart and tstop parameters are subselected, and other spin ephemerides are removed. See Appendix C for an explanation of validity of a spin ephemeris. Also, ephemeris remarks that are effective in (a part of) the time interval are subselected, and other ephemeris remarks are removed.
 If filter = SOLAREPH, the tool reads solareph parameter from the parameter interface of the tool, and filters the loaded pulsar ephemeris data by solar system ephemeris name. Spin ephemerides and orbital ephemerides that are derived using the solar system ephemeris specified by solareph parameter are subselected, and those using a different solar system ephemeris are removed.
 If filter = NONE, the tool does not filter the loaded pulsar ephemeris data.
 The tool reads author parameter from the parameter interface of the tool, in order to use its value in the next step for AUTHOR header keyword of the output pulsar ephemerides database.
 The tool creates an output pulsar ephemerides database file with a file name specified by outfile parameter, and writes all the remaining contents of the loaded pulsar ephemerides database into the output file.
Chapter 3. Details of common data processing procedures
Despite of differences in functionality and capability of the tools, most of them share common procedures such as pulsar ephemeris loading and photon arrival time corrections. This chapter describes common procedures of data processing by most of the pulsar analysis tools. Not all the steps are performed by all those tools, and not all the tools behave in the same way as others in all steps. Most of the tools, however, share procedures to prepare for their data processing, which this chapter describes in greater detail.
Each section of this chapter starts with a list of the tools that actually perform procedures described in the section. The tool names are listed after "Performed by:" tag at the beginning of each section. Also, it is summarized in the following tabulated matrix, where "Yes" means the step is performed by the tool, and "No" mean it is not. The matrix also outlines which part of the procedures each of the tools performs for your reference.
Step 1. Open event files
Performed by: gtpphase, gtophase, gtpsearch, gtpspec, gtptest
The tool evaluates evfile parameter, and creates a list of file names from the parameter value. If the value of evfile parameter starts with an atmark ("@"), the part that follows the atmark is interpreted as the name of a text file which contains a list of event file names, and the contents of the text file is set to the file name list. Otherwise, the value of evfile parameter is interpreted as an event file name, and set to the first and only element of the file name list.
Each file in the event file name list is then opened for an event extension and a Good Time Interval (GTI) extension in the order of appearance in the list. The name of an event extension to open is specified by evtable parameter, while a GTI extension name is fixed to "GTI".
When opening those extensions, the tool reads several FITS header keywords to determine a data format of each file. For example, the tool may read TELESCOP and INSTRUME header keywords to identify a mission name and an instrument type, as well as TIMEFRAM and TIMESYS header keywords to judge which arrival time corrections have already been applied to various kinds of times recorded in the file. If the data format is not supported by the tool, it produces an error message and terminates the process.
Step 2. Determine leap second file name
Performed by: gtbary, gtephem, gtpulsardb, gtpphase, gtophase, gtpsearch, gtpspec
The tool reads leapsecfile parameter from the parameter interface of the tool, and takes its value as the name of a leap second file to use in the later stage of the process, unless the value of leapsecfile parameter is DEFAULT. If leapsecfile parameter is set to DEFAULT, the tool looks for a default leap second file, although it is recommended to explicitly set leapsecfile parameter to the name of a leap second file of your choice if your analysis is sensitive to leap seconds, for the following reasons.
 An insertion or a removal of a leap second is unpredictable by nature, and a delivery of the latest leap second file does not necessarily synchronize with a software release. As a result, a default leap second file can be obsolete for your data set, and may not contain all the leap seconds necessary for your analysis.
 A default leap second file is a file named leapsec.fits under a directory specified by TIMING_DIR environment variable. However, using the environment variable to control a search for a leap second file could be misleading, because the value of the environment variable can be modified by a wrapper script, for instance, before the flow of execution reaches to the point where the tool reads it.
The tool needs a leap second file when UTC time system is involved in computation. Examples of such computations include, but are not limited to, interpreting a time expression in UTC and computing a time difference between two moments in time, one of which is expressed in UTC. The tool likely performs a UTCinvolved computation when timesys or usersys is set to UTC, or a pulsar ephemeris is given as defined in UTC.
Step 3. Interpret user's request on arrival time corrections
Performed by: gtpphase, gtophase, gtpsearch, gtpspec
The tool evaluates tcorrect parameter to interpret a user's request on which arrival time corrections should be applied to photon arrival times prior to computations of physical or statistical quantities from time series data, with an exception of gtophase which assumes that a user requests a predetermined combination of arrival time corrections.
There are three levels of request for each type of arrival time corrections: suppressed, allowed, or
required. If an arrival time correction is requested to be suppressed, the tool will not apply the arrival time correction on photon arrival times. If an arrival time correction is requested to be allowed, the tool searches for ephemeris data for the arrival time correction after loading all given pulsar ephemerides in the later stage of the process. If a set of ephemeris data is found for the arrival time correction, then the arrival time correction will be applied. Otherwise, it will not. If an arrival time correction is requested to be required, the tool searches for ephemeris data for the arrival time correction after loading all given pulsar ephemerides in the later stage of the process. If a set of ephemeris data is found for the arrival time correction, then the arrival time correction will be applied. Otherwise, the tool produces an error message and terminates the process.
The following table summarizes tool's interpretations of user's request on arrival time corrections for all possible values of tcorrect parameter. Due to differences in functionality and purposes of each tool, some tools have a slightly different way to interpret a user's request than others. Also in the table is gtophase, which does not have tcorrect parameter, showing that gtophase behaves as if a user requested arrival time corrections as indicated in the table.
Pulsar tool 
tcorrect value 
Barycentric correction 
Binary demodulation 
Pdot cancellation 
gtpphase 
NONE 
Suppressed 
Suppressed 
Suppressed 
AUTO 
Allowed 
Allowed 
Suppressed 
BARY 
Required 
Suppressed 
Suppressed 
BIN 
Required 
Required 
Suppressed 
ALL 
Required 
Required 
Suppressed 
gtophase 

Required 
Required 
Suppressed 
gtpsearch gtpspec 
NONE 
Suppressed 
Suppressed 
Suppressed 
AUTO 
Allowed 
Allowed 
Allowed 
BARY 
Required 
Suppressed 
Suppressed 
BIN 
Required 
Required 
Suppressed 
PDOT 
Required 
Suppressed 
Required 
ALL 
Required 
Required 
Required 
Note that not all possible combinations of arrival time corrections can be requested, because some combinations are physically inappropriate and should not be allowed. For example, because the barycentric correction must be applied before binary demodulation, a user cannot request to have the barycentric correction
suppressed and the binary demodulation required.
Step 4. Load pulsar ephemerides
Performed by: gtephem, gtpphase, gtophase, gtpsearch, gtpspec
The tool reads pulsar ephemerides and related remarks, and stores them in a memory space for later use. There are two types of pulsar ephemeris for the tool to load: a spin ephemeris and an orbital ephemeris. Naively speaking, a spin ephemeris is a set of parameters that describes temporal changes in pulsar's spin frequency, while an orbital ephemeris is a set of parameters that characterize propagation delays in a binary system that a pulsar of interest belongs to. Not all the pulsar analysis tools listed above loads the both types of pulsar ephemeris. An ephemeris remark can be any statement by a provider of pulsar ephemerides who wants to call user's attention, such as a report of a glitch or excessive timing noise during a certain period of time. For chatter = 4 or larger, if any information (spin ephemerides, orbital ephemerides, or ephemeris remarks) is loaded from pulsar ephemerides database file(s), the tool compiles a brief summary of the creation history of the loaded pulsar ephemerides database file(s), and reports the summary for human reading.
 Step 41. Load spin ephemerides

Performed by: gtephem, gtpphase, gtpsearch
The tool evaluates ephstyle parameter, and loads spin ephemerides as described in the following table. Because of the functionality and the capability of the tool, gtephem does not have ephstyle parameter, and loads spin ephemerides as if ephstyle were set to DB. For each spin ephemeris loaded, the tool chooses an appropriate numerical model of temporal changes in pulsar's spin frequency as described in Appendix A.
Pulsar tool 
ephstyle value 
Spin ephemerides are loaded... 
gtpphase gtpsearch 
DB 
from SPIN_PARAMETERS extension(s) of pulsar ephemerides database file(s) specified by psrdbfile parameter. Spin ephemerides are subselected by pulsar name given by psrname and, if requested, by solar system ephemeris name given by solareph parameter, before being stored in a memory space for later use. Subselection of spin ephemerides by solar system ephemeris name will be performed only when matchsolareph parameter is set to PSRDB or ALL. 
FREQ 
from parameters ephepoch, timeformat, timesys, ra, dec, phi0, f0, f1, and f2. For a tool without phi0 parameter in the parameter interface, phi0 parameter is assumed to be zero (0.0). 
PER 
from parameters ephepoch, timeformat, timesys, ra, dec, phi0, p0, p1, and p2. For a tool without phi0 parameter in the parameter interface, phi0 parameter is assumed to be zero (0.0). 
gtephem 

in the same way as gtpphase and gtpsearch do for ephstyle = DB. 
 Step 42. Load orbital ephemerides

Performed by: gtephem, gtpphase, gtophase, gtpsearch, gtpspec
Unless binary demodulation is requested to be suppressed (see Step 3 for types of request for arrival time corrections), orbital ephemerides are loaded from ORBITAL_PARAMETERS extension(s) of pulsar ephemerides database file(s) specified by psrdbfile parameter. Orbital ephemerides are subselected by pulsar name given by psrname and, if requested, by solar system ephemeris given by solareph parameter, before being stored in a memory space for later use. Subselection of orbital ephemerides by solar system ephemeris name will be performed only when matchsolareph parameter is set to PSRDB or ALL. For each orbital ephemeris loaded, the tool chooses an appropriate numerical model of pulsar's binary orbit as described in Appendix A.
 Step 43. Load ephemeris remarks

Performed by: gtephem, gtpphase, gtophase, gtpsearch, gtpspec
The tool evaluates reportephstatus parameter, and if the value of reportephstatus parameter is a logical true, then the tool loads ephemeris remarks from REMARKS extension of pulsar ephemerides database file(s) specified by psrdbfile parameter.
Step 5. Initialize arrival time corrections
Performed by: gtpphase, gtophase, gtpsearch, gtpspec
The tool determines which arrival time corrections to apply to photon arrival times on the fly, defines a time axis for time series analysis, and initializes arrival time corrections that are determined to be applied. If chatter parameter is set to 3 or larger, the tool summarizes its determinations on arrival time corrections, and reports the summary for human reading.
 Step 51. Determine arrival time corrections to apply

Performed by: gtpphase, gtophase, gtpsearch, gtpspec
The tool determines which arrival time corrections to apply to photon arrival times by a combination of a user's request and availability of ephemeris data for requested arrival time corrections. The following table summarizes tool's determination on whether each type of arrival time correction should be applied.
Arrival time correction 
Pulsar tool 
User's request 
The arrival time correction is determined... 
binary demodulation 
gtpphase gtophase gtpsearch gtpspec 
suppressed 
not to be applied. 
allowed 
to be applied if and only if at least one orbital ephemeris is loaded. 
required 
to be applied if at least one orbital ephemeris is loaded. If none is loaded, the tool produces an error message and terminates the process. 
pdot cancellation 
gtpsearch 
suppressed 
not to be applied. 
allowed 
to be applied if and only if at least one spin ephemeris is loaded. 
required 
to be applied if at least one spin ephemeris is loaded. If none is loaded, the tool produces an error message and terminates the process. 
gtpspec 
suppressed 
not to be applied. 
allowed required 
to be applied. 
barycentric correction 
gtpphase gtophase gtpsearch gtpspec 
suppressed 
not to be applied. If binary demodulation or pdot cancellation are determined to be applied, the tool produces an error message and terminates the process. 
allowed required 
to be applied. 
 Step 52. Define time axis for time series analysis

Performed by: gtpsearch, gtpspec
The tool determines a time origin and a time system for time series analysis and stores in a memory space for later use. There are multiple ways to specify a time origin, and a user can choose a desired one by setting an appropriate value to timeorigin parameter.
timeorigin value 
A time origin is set to... 
START 
the start time of a given data set. See Appendix B for definition of the start time of a data set. 
STOP 
the stop time of a given data set. See Appendix B for definition of the stop time of a data set. 
MIDDLE 
the midpoint of the start time and the stop time of a given data set. See Appendix B for definition of a start and a stop time of a data set. The midpoint of the two moments in time is computed in a time system specified by TIMESYS header keyword of the given event file, or the first event file in the event file list given through @mark convention if more than one event files are given. 
USER 
a moment in time that usertime parameter points to as in a time format specified by userformat parameter in a time system specified by usersys
parameter. 
A time system for time series analysis is determined depending on arrival time corrections (barycentric correction, binary demodulation, and pdot cancellation) that are requested. If at least one of them is requested, TDB system is chosen for time series analysis. If no arrival time corrections are requested, the tool reads a TIMESYS header keyword value from all the given event files, and checks whether they are all identical to each other. If they are identical to each other, then the time system is chosen for time series analysis. Otherwise, the tool produces an error message and terminates the process.
 Step 53. Initialize pdot cancellation

Performed by: gtpsearch, gtpspec
If pdot cancellation is determined to be applied in the earlier section of this step, (see above for the determination procedure), the tool determines values of the following parameters for pdot cancellation and a time system in which pdot cancellation computations are performed, and stores them in a memory space for later use. If pdot cancellation is determined not to be applied, the tool does nothing in this section.
r_{1}: ratio of the first timederivative of spin frequency over the spin frequency at the time origin
r_{2}: ratio of the second timederivative of spin frequency over the spin frequency at the time origin
Due to differences in functionality and capability of the tools, each
tool has its own way to determine them.
gtpsearch chooses one spin ephemeris for the time origin from the loaded spin ephemerides as described in Appendix C. Validity of a chosen spin ephemeris is not required. See Appendix C for an explanation of validity of a spin ephemeris. Using the chosen spin ephemeris, the tool computes spin frequency and its first and second timederivatives at the time origin (see Appendix A for how to compute those quantities). The tool then computes a ratio of the first derivative over the spin frequency and another ratio of the second derivative over the spin frequency, and takes them as r_{1} and r_{2}. Finally, the tool identifies an appropriate time system in which ephemeris computations with the chosen spin ephemeris are performed (see Appendix A for how to determine a time system for ephemeris computations), and takes it as a time system to be used for pdot cancellation.
 gtpspec evaluates ephstyle parameter. If ephstyle = FREQ, the tool reads f1f0ratio and f2f0ratio parameters from the parameter interface of the tool, and takes them as r_{1} and r_{2}, respectively. In other words, r_{1} and r_{2} are set as:
r_{1} = f1f0ratio
r_{2} = f2f0ratio
where f1f0ratio and f2f0ratio are parameter values taken from the parameter interface of the tool. If ephstyle = PER, the tool reads p1p0ratio and p2p0ratio parameters from the parameter interface of the tool, and computes r_{1} and r_{2} as:
r_{1} =  p1p0ratio
r_{2} = 2 × (p1p0ratio)^{2}  p2p0ratio
where p1p0ratio and p2p0ratio are parameter values taken from the parameter interface of the tool. In both cases, the tool performs pdot cancellation computations in a time system for time series analysis determined in Step 52.
 Step 54. Initialize barycentric correction

Performed by: gtpphase, gtophase, gtpsearch, gtpspec
The tool reads scfile, sctable, and
angtol parameters from the parameter interface of the tool, and stores them in a memory space for later use. See Appendix D for how the tool uses those parameter values. The tool also evaluates matchsolareph parameter, and if matchsolareph = EVENT or ALL, then the tool checks all event files in which event times are already barycentric, in order to enforce consistent use of a single solar system ephemeris in barycentric corrections throughout all the given event file(s). For other values of matchsolareph parameter, the tool does not perform the checking. For each event file whose event times are already barycentric, the tool performs the checking as follows.
 The tool identifies the name of solar system ephemeris with which barycentric times in the event file were computed. Typically, it is recorded as a PLEPHEM header keyword value in the event file.
 The tool compares it with the value of solareph parameter. If they match, the tool continues processing the event file. Otherwise, the tool produces an error message and terminates the process.
In addition, gtophase and gtpspec reads ra and
dec parameters from the parameter interface of the tool, and stores them in a memory space for later use as a direction in the sky for which barycentric corrections will be computed.
 Step 55. Initialize binary demodulation

Performed by: none
The tool does nothing to initialize binary demodulation.
Step 6. Summarize ephemeris status
Performed by: gtpphase, gtophase, gtpsearch, gtpspec
The tool examines the ephemeris information that are loaded in Step 4 (i.e., spin ephemerides, orbital ephemerides, and ephemeris remarks), and produces a warning message if any noteworthy ephemeris status is identified. A warning message produced in this step is only for human reading, and does not affect data processing by the tool at all. It is designed merely to notify a user availability of ephemeris data throughout anticipated data processing (i.e., continuity of ephemeris coverage without a gap), or to call user's attention to a certain special condition such as a glitch.
The loaded ephemeris information is examined within a time interval derived from a time coverage of given event data. First, the tool computes the start time of a given data set. See Appendix B for definition of the start time of a data set. Then, the start time is applied arrival time corrections that are determined to be applied to photon arrival times in Step 51. See Appendix D for how the tool applies arrival time corrections. If it is successful, then the tool starts ephemeris examination at the corrected start time. If it is unsuccessful (which can happen due to lack of ephemeris data for the moment in time, for example), the tool starts ephemeris examination at the uncorrected start time. Similarly, the tool determines a stopping time of ephemeris examination based on the stop time of a given data set.
Note that unsuccessful arrival time corrections in this step do not produce an error message or a warning. That is because unsuccessful arrival time corrections on the start time and the stop time do not necessarily mean arrival time corrections will be unsuccessful in the later process of data.
Also note a warning message from the ephemeris examination performed in this step may not be complete. If arrival time corrections are unsuccessful in computing the start time or the stop time of the examination time interval, ephemeris gap(s) or ephemeris remark(s) may be missed near the uncorrected time(s), for example.
The tool then examines the loaded ephemeris information within the time interval determined above, and reports the result of examination. Not all the tools examine all kinds of ephemeris status. Depending on the tool's functionality, each tool has a set of ephemeris status to examine, as described below.
 gtpphase examines and reports ephemeris gaps within the time interval and ephemeris remarks that overlaps with the time interval.
 gtophase, gtpsearch, and gtpspec examines and reports only ephemeris remarks that overlaps with the time interval.
For each ephemeris status reported in a warning message, the tool also reports a time interval during which the ephemeris status persists. Based on a warning message produced in this step, a user may choose to further subselect events by event time, in order to avoid potential problem(s) in user's analysis caused by the reported ephemeris status.
Appendix
Appendix A. Numerical models of pulsar's timing properties
Pulsar's timing properties reflect two types of pulsar's physical motions: a rotational motion and an orbital motion. In a conventional time series analysis, a rotational motion of a pulsar is characterized by a numerical model that describes temporal changes in pulsar's spin frequency. There are a few models for temporal changes in spin frequency. Each spinfrequency model is parametrized by a set of parameters in its own way, and a spin frequency and its time derivatives at an arbitrary time can be computed with the parameter set. Also, it is common in a conventional time series analysis to characterize a pulsar's orbital motion by a numerical model that describes a propagation time for a photon that is emitted at a pulsar and travels in a binary system that the pulsar belongs to. Various physical and empirical models have been developed for an orbital motion of a pulsar in a binary system. Each binary model is parametrized by a set of parameters in its own way, and a photon propagation time in a binary system can be computed with the parameter set.
This section lists the spinfrequency models and the binary orbital models that are implemented in the pulsar analysis tools, and explains how the pulsar analysis tools compute relevant quantities for time series analysis with given spin and orbital parameters.
 Ephemeris computations with spin parameters

A spinfrequency model can be represented by a function of time that computes the number of revolutions that the pulsar will turn until a given point of time since a reference time. When the tool loads spin ephemerides, the tool identifies an appropriate model of spinfrequency history for each spin ephemeris, and chooses a function of time to compute a pulsar position, a pulse phase, a pulse frequency, and time derivatives of pulse frequency whenever necessary in data processing.
Letting α(t) and δ(t) be the pulsar's Right Ascension and Declination, respectively, N(t) the number of pulsar's revolutions at time t, φ(t) a pulse phase at time t, f(t) a pulse frequency at time t, and f^{(n)}(t) the nth time derivative of pulse frequency at time t, φ(t) is computed as the fractional part of N(t) independent of spinfrequency model. By convention, φ(t) is computed to be greater than or equal to zero (0.0) and less than unity (1.0). Other quantities are computed as described below, depending on a spinfrequency model that the tool chooses for a given spin ephemeris.
 If a spin ephemeris is loaded from a SPIN_PARAMETERS extension of a pulsar ephemerides database file with header keyword EPHSTYLE = "FREQ", the quantities are computed in TDB time system as:
α(t) = RA
δ(t) = DEC
N(t) = φ_{0} + F0 × (t  t_{0}) + 1/2 × F1 × (t  t_{0})^{2} + 1/6 × F2 × (t  t_{0})^{3}
f(t) = F0 + F1 × (t  t_{0}) + 1/2 × F2 × (t  t_{0})^{2}
f^{(1)}(t) = F1 + F2 × (t  t_{0})
f^{(2)}(t) = F2
f^{(n)}(t) = 0 for n ≥ 3
where RA, DEC, F0, F1, and F2 are field values of FITS columns of such names, φ_{0} a pulse phase constant which is further explained below, t_{0} a moment in time that the sum of EPOCH_INT and EPOCH_FRAC points to as an MJD number in TDB time system. Pulse phase constant φ_{0} is the fractional part of N_{0}, which is then computed as
N_{0} =  F0 × (t_{a}  t_{0})  1/2 × F1 × (t_{a}  t_{0})^{2}  1/6 × F2 × (t_{a}  t_{0})^{3},
where t_{a} is a moment in time that the sum of TOABARY_INT and TOABARY_FRAC points to as an MJD number in TDB time system. Combining the definitions above, one obtains φ(t_{a}) = 0.0 (zero). In other words, a pulse phase at a time of pulse arrival is always zero (0.0), allowing multiple sets of spin parameters to be used to assign pulse phase values to photon arrival times in a consistent manner, to the extent of reliability of ephemeris data used.
 If a spin ephemeris is loaded from the parameter interface of the tool with ephstyle = FREQ, the quantities are computed in a time system specified by timesys parameter as:
α(t) = ra
δ(t) = dec
N(t) = phi0 + f0 × (t  t_{0}) + 1/2 × f1 × (t  t_{0})^{2} + 1/6 × f2 × (t  t_{0})^{3}
f(t) = f0 + f1 × (t  t_{0}) + 1/2 × f2 × (t  t_{0})^{2}
f^{(1)}(t) = f1 + f2 × (t  t_{0})
f^{(2)}(t) = f2
f^{(n)}(t) = 0 for n ≥ 3
where ra, dec, phi0, f0, f1, and f2 are parameter values taken from the parameter interface of the tool, t_{0} a moment in time that ephepoch parameter points to as in a time format specified by timeformat parameter in a time system specified by timesys parameter.
 If a spin ephemeris is loaded from the parameter interface of the tool with ephstyle = PER, the quantities are computed in a time system specified by timesys parameter as:
α(t) = ra
δ(t) = dec
N(t) = phi0 + I(t  t_{0}, p0, p1, p2)
f(t) = 1 / q_{0}
f^{(1)}(t) =  q_{1} / q_{0}^{2}
f^{(2)}(t) =  q_{2} / q_{0}^{2} + 2 q_{1}^{2} / q_{0}^{3}
f^{(n)}(t) = D_{n}(q_{0}, q_{1}, p2)
where ra, dec, phi0, p0, p1, and p2 are parameter values taken from the parameter interface of the tool, t_{0} a moment in time that ephepoch parameter points to as in a time format specified by timeformat parameter in a time system specified by timesys parameter, I(t, t_{0}, p0, p1, p2) the definite integral of f(t) over time range from t_{0} to t, q_{0} and q_{1} are a pulse period and its first time derivative at time t as defined by:
q_{0} = p0 + p1 × (t  t_{0}) + 1/2 × p2 × (t  t_{0})^{2}
q_{1} = p1 + p2 × (t  t_{0})
and finally D_{n}(q_{0}, q_{1}, p2) the nth time derivative of f(t) expressed as a polynomial of q_{0}, q_{1}, and p2. The explicit forms of the functions I(t, t_{0}, p0, p1, p2) and D_{n}(q_{0}, q_{1}, p2) are derived and described in a separate document titled "Computing pulse phase and pulse frequency using pulse period ephemeris" (PDF file, 92 kB) for your reference.
 Ephemeris computations with orbital parameters

An orbital model of a binary system can be represented by a function of time that computes a time delay in a photon arrival time caused by an orbital motion of a pulsar in a binary system. When the tool loads orbital ephemerides, the tool identifies an appropriate binary orbital model for each orbital ephemeris, and chooses a function of time to compute a photon propagation time in a binary system. Based on the chosen function, the tool computes a photon arrival time at which a photon of interest would have arrived at the spacecraft if the photon left at the center of gravity of the binary system, instead of the pulsar. This computation, called binary demodulation, is performed as a part of photon arrival time corrections, if so requested. Binary demodulation is an iterative process to search for a corrected photon arrival time (see Appendix D for details), and each orbital model defines the maximum tolerance in time difference as convergence criteria in the iterative search.
 If an orbital ephemeris is loaded from a ORBITAL_PARAMETERS extension of a pulsar ephemerides database file with header keyword EPHSTYLE = "DD", time delays are computed in TDB time system based on the DamourDeruelle model (Damour and Deruelle 1985 and 1986) with no aberration (i.e., A = B = 0). The explicit formula of time delay was originally derived by Damour and Deruelle (1986), and can be found as Equation 7 by Taylor and Weisberg (1989), for example, where a brief description of the model and an excellent summary of other binary models are also given. The maximum tolerance in time difference for binarydemodulation is set to 1 nanosecond (1 × 10^{9} seconds).
References
Damour, T., and Deruelle, N. 1985, Ann. Inst. H. Poincare (Physique Theorique), 43, 107.
Damour, T., and Deruelle, N. 1986, Ann. Inst. H. Poincare (Physique Theorique), 44, 263.
Taylor, J. H., and Weisberg, J. M. 1989, ApJ 345, 454.
Appendix B. Time boundaries of event data
Time boundaries of a give set of data often plays an important roll in time series analysis. In the pulsar analysis tools, the time boundaries of given event data are defined as described below.
 Start time
 a moment in time that a START column value of the first FITS row in GTI extension of the given event file points to. If more than one event files are given, the tool computes the start time for each event file, and takes the earliest one among them as the start time for the entire data set.
 Stop time
 a moment in time that a STOP column value of the last FITS row in GTI extension of the given event file points to. If more than one event files are given, the tool computes the stop time for each event file, and takes the latest one among them as the stop time for the entire data set.
Appendix C. Best ephemeris for a given moment in time
In several occasions in data processing by the pulsar analysis tools, the tool needs to choose a spin ephemeris or an orbital ephemeris that is best used to compute a physical quantity (such as a spin frequency or a photon propagation delay in a binary system) at a given moment in time. Which ephemeris is best among the loaded ephemerides depends on a quantity to compute and a purpose of the computation, as well as a moment in time for which an ephemeris is chosen. This section describes the ephemeris choosing procedures that are implemented in the pulsar analysis tools.
 Best spin ephemeris when validity required

Each spin ephemeris comes with a time interval during which temporal changes of spin frequency is well reproduced by the ephemeris. The time interval is called a validity window of the spin ephemeris, and the spin ephemeris is considered valid during the time interval. When validity of a chosen ephemeris is required, the tool chooses a spin ephemeris from the loaded spin ephemerides in the following procedure.
 The tool searches for a spin ephemeris which is valid at the given moment in time.
 If only one spin ephemeris is found, the tool chooses it as the best spin ephemeris.
 If more than one spin ephemerides are found, the tool takes the following steps to compare them until only one spin ephemeris is left. Then, the tool chooses it for the best spin ephemeris.
 Ephemeris with a validity window that starts latest
 Ephemeris with the longest validity window
 Ephemeris that appears latest in the loaded spin ephemerides
 If none is found, the tool reports an error.
 Best spin ephemeris when validity not required

When validity of a chosen ephemeris is not required, the tool chooses a spin ephemeris from the loaded spin ephemerides in the following procedure.
 The tool tries to find the best spin ephemeris which is valid at the given moment in time in the procedure described above.
 If found, the tool chooses it as the best spin ephemeris.
 If none is found (i.e., the tool reports an error in finding the best spin ephemeris which is valid at the given moment in time), the tool searches for a spin ephemeris that is closest to the given moment in time in the following procedure.
 For each of the loaded spin ephemeris, the tool computes in TDB system the time difference between the start time of its validity window and the given moment in time, and the time difference between the stop time of its validity window and the given moment in time.
 The tool compares the lengths of the two time differences, and takes the shorter one as the time distance between the spin ephemeris and the given moment in time.
 The tool searches for a spin ephemeris with the shortest time distance to the given moment in time, and takes it as the closest spin ephemeris to the given moment in time. If more than one spin ephemerides are tied, the tool takes the one that appears latest in the loaded spin ephemerides.
 The tool chooses the closest spin ephemeris to the given moment in time for the best spin ephemeris.
 Best orbital ephemeris

Unlike a spin ephemeris, an orbital ephemeris does not come with a validity window, based on which the best spin ephemeris can be defined for an arbitrary moment in time. Instead, the tool defines a reference epoch for each orbital model, and uses it to choose the best orbital ephemeris among the loaded orbital ephemerides as described below. Unless otherwise noted, the tool uses a time of periastron passage as a reference epoch of an orbital ephemeris, since it is commonly seen in orbital parameters of most binary models.
 For each of the loaded orbital ephemeris, the tool computes in TDB system the time difference between the reference time of the orbital ephemeris and the given moment in time, and takes it as the time distance between the orbital ephemeris and the given moment in time.
 The tool searches for an orbital ephemeris with the shortest time distance to the given moment in time, and takes it as the closest orbital ephemeris to the given moment in time. If more than one orbital ephemerides are tied, the tool takes the one that appears latest in the loaded orbital ephemerides.
 The tool chooses the closest orbital ephemeris to the given moment in time for the best orbital ephemeris.
Appendix D. Arrival time corrections
Three arrival time corrections are implemented for the pulsar analysis tools: the barycentric correction, binary demodulation, and pdot cancellation. When requested, the tool applies those corrections to a given moment in the order listed above. This section describes each of the arrival time corrections in detail.
 Barycentric correction

For an event time (or any other representation of time) in an event file, the tool computes a barycentric time, which is a moment in time when the photon would have arrived at the solar system barycenter, as described below.
 The tool computes the pulsar position (Right Ascension and Declination) in the sky for a given moment in time. gtpphase and gtpsearch chooses the best spin ephemeris for the given moment in time (see Appendix C for ephemeris choosing procedures), and computes the pulsar position at the given moment in time with the spin ephemeris (see Appendix C for ephemeris computation procedures). gtophase, gtpspec, and gtbary read ra and dec parameters in the parameter interface of the tool, and use them as the sky position for which barycentric corrections are applied.
 If it is indicated that the event time is already a barycentric time (typically by TIMEREF header keyword of the event file being SOLARSYSTEM), the tool checks consistency in pulsar position, by comparing the pulsar position for which the barycentric time was computed with the one for which the barycentric corrections would have been applied to the event time if it were not barycentric. If a separation angle between the two pulsar positions is less than or equal to the value of angtol parameter, then the tool determines they coincide with each other. Otherwise, the tool produces an error message and terminates the process.
 The tool reads spacecraft data from a table (usually a FITS extension) specified by sctable parameter in a spacecraft file specified by scfile parameter, and computes the spacecraft position in reference to the center of the Earth at the given moment in time.
 The tool obtains or computes the following quantities at the given moment in time.
 r_{obs} : the vector pointing from the solar system barycenter to the spacecraft
 n_{psr} : a unit vector pointing to the pulsar
 T_{sys} : the time difference between TT and TDB systems (TDB  TT).
 r_{sat} : the vector pointing from the geocenter (the center of the Earth) to the spacecraft
 v_{earth} : the velocity of the Earth with respect to the solar system barycenter
 θ : the pulsarSunEarth angle, i.e., the angle formed between a line connecting the pulsar and the center of the Sun and a line connecting the Earth and the center of the Sun. Note it is the center of the Sun to use here, not the solar system barycenter.
 The tool computes a barycentric time for the given moment in time. Let t_{obs} be the given moment in time, t_{bary} a barycentric time for it, then the tool computes it as:
t_{bary} = t_{obs} + Δt_{g} + Δt_{E}  Δt_{S}
where Δt_{g} is a geometric time delay in the solar system (i.e., a light travel time without gravitational effects), Δt_{E} is the Einstein corrections, Δt_{S} is the Shapiro delay due to the gravitational field of the solar system. Each term on the righthand side is computed as:
Δt_{g} = r_{obs} · n_{psr} / c
Δt_{E} = T_{sys} + r_{sat} · v_{earth} / c^{2}
Δt_{S} =  (2GM_{sun} / c^{3}) log(1 + cosθ)
where c is the speed of light, G the gravitational constant, and M_{sun} the solar mass.
 Binary demodulation

For a barycentric time, the tool computes a binarydemodulated time, which is a moment in time when the photon would have arrived at the solar system barycenter if it were emitted at the center of gravity of the binary system that the pulsar belongs to. The tool performs binary demodulation on the fly if it is determined to be applied in Step 51 of Chapter 3. Binary demodulation is an arrival time correction for an orbital delay, i.e., a time delay caused by the fact that the photon was emitted and traveled within a binary system, with respect to the case where the photon were emitted at the center of the gravity of the binary system. Note that computations of orbital delays are dependent on a binary orbital model, and may or may not include gravitational effects depending on a binary orbital model.
As described below, the procedure is iterative, but typically a convergence is achieved in a few iterations.
 The tool computes initial guesses on a binarydemodulated time and an orbital delay for the photon.
 The tool takes the given barycentric time as an initial guess on a binarydemodulated time.
 The tool chooses the best orbital ephemeris for the given barycentric time.
 With the chosen orbital ephemeris, the tool computes an orbital delay for a photon emitted at the given barycentric time, and takes it as an initial guess on an orbital delay.
 The tool iterates the following procedures until it converges (the convergence criteria is also described below). In the first iteration, the tool takes the initial guess of a binarydemodulated time as a candidate time, and the initial guess of an orbital delay as a tentative orbital delay.
 The tool subtracts the tentative orbital delay from the candidate time, and takes the result as a new candidate time.
 The tool chooses the best orbital ephemeris for the new candidate time.
 With the chosen orbital ephemeris, the tool computes an orbital delay for a photon emitted at the new candidate time, and takes it as a new tentative delay.
 The tool adds the new tentative delay to the new candidate time, and compares the result with the given barycentric time. If the time difference is less than or equal to the maximum tolerance in time difference for the chosen binary orbital model (see Appendix A for details), it is deemed converged. Otherwise, take another iteration with the new candidate time and the new tentative orbital delay.
 The tool takes the candidate time at the convergence as a binarydemodulated time.
 Pdot cancellation

The tool modifies a barycentric time or a binarydemodulated time to cancel or reduce the effect of frequency derivatives by the given amount, in order to make potential pulsations in the time series data appear to have a constant frequency over time.
Letting t_{unc} be a photon arrival time that is to be corrected for pdot cancellation, t_{pdot} a corrected time, then the tool computes it in a time system determined in Step 53 of Chapter 3 as:
t_{pdot} = t_{unc} + 1/2 × r_{1} × (t_{unc}  t_{origin})^{2} + 1/6 × r_{2} × (t_{unc}  t_{origin})^{3}
where r_{1} and r_{2} are the pdot cancellation parameters as defined in Step 53 of Chapter 3, and t_{origin} is the time origin that is computed in Step 52 of Chapter 3.