The Planetary Ephemeris Program: Capability, Comparison, and Open Source Availability

This paper provides the first published description of the Planetary Ephemeris Program (PEP), which has been in use and evolving for nearly six decades. PEP is a general-purpose astrometric data-analysis program that can compute numerical ephemerides and can analyze simultaneously a heterogeneous collection of astrometric data. Very few programs like PEP exist and, to our knowledge, PEP is the only one that has been freely available, including source code, to interested parties. Among the other programs with comparable scope are DPODP (JPL, USA; Park et al. 2021), INPOP (IMCCE Obs. Paris, France; Viswanathan et al. 2018), the analysis model at the Institut für Erdmessung, Leibniz Universität Hannover, Germany (Müller et al. 2019), and EPM (Laboratory of Ephemeris Astronomy, Russia; Pitjeva & Pitjev 2014).

Ephemeris and modeling codes such as PEP are flexible tools that have played a critical role in determining physical parameters of solar system bodies such as orbital initial conditions, planetary masses, and the length scale (the au) in terrestrial length units (Ash et al. 1967), and in testing fundamental physics theories. PEP has a particularly rich history with the latter. It has been used to measure general relativistic effects such as the de Sitter precession of the moon (Shapiro et al. 1988) and the Shapiro time delay (Reasenberg et al. 1979). In addition, PEP has been used to test physics theories that go "beyond the standard model" (Battat et al. 2007, 2008).

Over the past fifty years, the suite of solar system astrometric data sets has grown in scope and precision, and advancements from the theory community have triggered an evolution in the types of scientific questions that we ask of the data. In response to these changes, PEP has been updated, both in terms of the types of observables it can analyze and in terms of the precision with which it can analyze observables. A major current effort with PEP is the analysis of the Lunar Laser Ranging (LLR) data from the Apache Point Observatory Lunar Laser-ranging Operation (APOLLO; Murphy et al. 2008, 2012) experiment in conjunction with the other sets of data described in Table 1. APOLLO represents an order-of-magnitude improvement in the precision of LLR, yielding an Earth–Moon distance equivalent of 1.2 mm median uncertainty in a single one hour observing session (Battat et al. 2009; Liang et al. 2017). As is frequently true of the analysis of new (improved) astrometric data, the model of the motion and rotation of the Moon needed to reduce these data to near their noise level likely has more free parameters than can be estimated simultaneously. This is in general true of the solar system model and data set. Under these conditions, it may be useful to add auxiliary data types that likely will reduce the degeneracy of the estimator and thus allow for the use of a more complex model for the rotation and center-of-mass motion of the Moon. Therefore, driven by this advance in LLR, PEP has been upgraded to handle pulse time-of-arrival (TOA) data from more types of pulsars than before, such as those in binary systems.

Such an upgrade benefits other analyses too and opens the possibility of new ones. For example, pulsar TOA data can be used in conjunction with VLBI data to tie the PEP's planetary reference frame to the extra-galactic radio reference frame (Bartel et al. 1996). Moreover, even without a frame tie, the observation of pulsars well away from the ecliptic offers a breaking of the degeneracy implicit in relying on data restricted to the vicinity of a single plane. For another example, the Dvali–Gabadadze–Porrati (DGP) model of gravity (Dvali et al. 2000), which was introduced to explain cosmic acceleration without dark energy, predicts a common-mode precession of all orbiting bodies in the solar system by a fixed angular rate. For coplanar orbits, this looks like a net rotation. The inclinations of the planetary orbits relative to the ecliptic break this degeneracy, but the small inclination angles strongly suppress the sensitivity to the predicted effect in DGP theory (Battat et al. 2008). The out-of-plane perspective provided by pulsar timing would enhance the sensitivity to DGP-like anomalies.

In Section 2, we provide an overview of PEP, including a description of the astrometric data sets it can process. This section concludes with a discussion of the PEP's availability, both historically and for the foreseeable future. In Section 3, we discuss the handling of pulsar timing observables in PEP and present the results of a validation of this capability against Tempo2 (Edwards et al. 2006; Hobbs et al. 2006), a program specialized for the analysis of TOA observations of pulses transmitted from a pulsar and in wide use by the pulsar community. Such comparisons of scientific codes are an essential part of the effort to ensure the validity of data analysis in an era with very large precision codes and evolving demands on these codes. A brief summary of the key points of the paper and some further observations are contained in Section 4. Finally, in the appendix, we provide the algorithm that allows us to make numerous numerical solution estimates at the end of a data analysis (when the estimate is converged and most nearly linear). Unlike most of the algorithms embedded in PEP, we believe this one does not have a counterpart in other analysis packages.

The PEP was originally created at MIT Lincoln Laboratory starting in the early 1960s. It was brought to the main campus half a decade later, which then evolved to be the center of PEP activity. In the early 1980s that center moved to the (then named) Harvard-Smithsonian Center for Astrophysics (CfA). The PEP is a general-purpose astrometric data-analysis program and can analyze simultaneously various kinds of astrometric data, as described in Table 2. That table includes the highest precision of the data analyzed for each category, which sets the level of model complexity (completeness) required. A top-level outline of PEP is provided in Table 3.

As might be expected from its inception era, the PEP is written in Fortran, albeit modernized to make use of language features not available to the original Fortran II implementation. In the current era, we use the GNU G Fortran compiler. In early versions of PEP, there were some calculations performed in "extended precision" using calls to assembly language routines. With GNU G Fortran, many of the calculations are done using IEEE-format 80-bit floating point numbers. Recently, we have been looking at other extended precision options again for some sensitive calculations. Over the decades, PEP has accreted new capabilities as needed for its planned uses. In Section 3, we discuss one such addition, the inclusion of the analysis of the time-of-arrival observations of pulses transmitted from a pulsar. ⁵ Although this addition uses algorithms generated within the PEP group, we are verifying that these algorithms are consistent with the best models used in the field. To that end, we have compared the PEP analysis of a particular set of pulsar data with the corresponding analysis made with Tempo2.

Over the past half century or so, the PEP has been used in over 20 doctoral theses including, for example, several on the fourth test of general relativity in the late 1960s and 1970s and in about 90 scientific papers including about 30 on tests of general relativity, competing theories and related cosmological issues such as a possible secular variation of the Newtonian "constant" of gravitation. These range from tests of the Shapiro time delay during the early years (Shapiro et al. 1968, 1971) to more recent papers on constraints on the DGP braneworld theory of gravity (Battat et al. 2008) and constraints on Standard-Model Extension Parameters (Battat et al. 2007).

The maintainers of the PEP have always made the PEP code available to members of the scientific community, thus starting well before the popular use of the term "open source." Those interested generally had an affiliation with the then current users and were directed to a member of the PEP community who would provide a copy of the code and some guidance in its use. We have since shifted to a more accessible and user-driven form of distribution. The PEP source code, utility programs, documentation, and "Unit Tests" are now all publicly available via GitHub, and are distributed under a Creative Commons Attribution-NonCommercial-ShareAlike license, making PEP open source in the modern sense (Battat 2021b). The Unit Tests provide examples of PEP run streams, which can conveniently be used as starting points for the development of new run streams, see Table 4. The table lists the operations performed in each of the tests. The collection of tests has grown along with PEP over the decades as new features were added to the latter. There are many technical memoranda written at MIT and CfA that describe algorithms and other key aspects of PEP. Several of these are available via GitHub and the digital collection continues to grow (Battat 2021a). Historically, PEP ran under various versions of IBM System 360 and System 370, then migrated to Sun OS, which became Solaris, then jumped to Linux and thus to basically any type of Unix, including macOS. It has been run under Windows experimentally, but not for production, as far as we know.

As for the pulsar TOA measurements, we also use NPs, described below, each derived from observation of a pulsar, usually over the order of an hour or less, rather than the raw data. In creating these NPs, the received signal from the pulsar is recorded at a variety of frequencies. In the Parkes Pulsar Timing Array Program, observations are conducted in the 20 cm (1400 MHz) band, or in the 10 cm (3100 MHz) and 40 cm (700 MHz) bands simultaneously (Manchester et al. 2013). Often, the signal is recorded simultaneously in two bands, either widely spaced (e.g., two selected from these three major bands) or closely spaced (e.g., two narrow bands both within the 20 cm band and separated only by about 60 MHz to avoid a band of interference). The number of observing sessions, each about 10–60 minutes long, may range from one to ten in a day and the typical observing cadence is two weeks. To generate an NP, the individual pulse waveforms are combined via integrating over time and frequency within one of the bands used during the observing session. The stability of the pulsar allows such integrations to be made easily. The resulting combined (averaged or folded) waveform is compared to a template, developed from a much larger ensemble of waveforms from the same pulsar and observing band, to determine the best estimate of the arrival time of a single representative pulse near the middle of an observing session and the middle of the frequency band over which the integration extended. This is the NP for that band and observing session.

We have recently upgraded the PEP code for the TOA of the radio (or, in principle, other electromagnetic emission) pulses from a pulsar. The formulation of this observable has been published for the pulsar analysis packages Tempo2 and, more recently, PINT (Luo et al. 2021). However, for our demonstration analysis, we chose to make only minor modifications to our decades-old simpler model that would work well enough for this purpose. We compared PEP to Tempo2 for real data, in particular, a set of 1065 NPs from PSR J1909-3744, spanning about 8 years (from 2002 December to 2011 February). ⁶ There was no consultation or coordination between the PEP and Tempo2 teams on algorithms used, except for the plasma, the reduction of station timing to UTC, and the treatment of the pulsar companion. PSR J1909-3744, which is part of the binary system, was discovered in late 2001 (Jacoby et al. 2003) at the 64-meter-diameter Parkes radio telescope in Australia (see Manchester et al. 2013). The inter-pulse period is 2.95 ms. The pulse itself has a full width at half maximum of 43 μs, with a sharp peak, and is extremely stable. These properties make this pulsar an excellent candidate for a software comparison.

3.1. Analysis of Data

3.2. Comparison of Analyses

As can be seen in Table 5, the two programs produce nearly the same formal standard deviations for all 10 parameters to the number of places shown. In three cases (dispersion rate, pulse offset, and pulse period rate), they differ by one in the least significant digit, which we do not take to be cause for concern. Similarly, the two programs produce nearly the same parameter estimates. In all but one case, the absolute difference divided by the standard deviation is under or well under 1. The exception is the pulse period, which is discussed below.

It is worth noting that the parameter estimates are all sensitive to context. For example, the adoption of a particular planetary ephemeris establishes the reference frame in which the position, proper motion, and parallax of the pulsar are measured. Accordingly, it was essential for us to use the same ephemeris in order to obtain comparable results. Needless to say, it was also obligatory to adopt matching reference epochs for the time-varying quantities (position, dispersion, and pulse phase). Indeed, we used a single reference epoch for all three. One remaining complication is the timescale. Tempo2 calculations are expressed in terms of the TCB (Barycentric Coordinate Time) scale, while PEP and the DE421 ephemeris use TDB (Barycentric Dynamical Time), and the two scales progress at different rates. The pulse period is the only parameter among the ten we compared that is directly affected by this complication. To make the comparison in Table 5, we converted the Tempo2 result for the pulse frequency expressed in pulses per TCB second into a pulse period in TDB seconds, using a scale factor (L_B ) extracted from Tempo2; but there are several other conventional values for L_B with significant variations relative to the precision required for our comparison. Since this parameter is the only one with an apparently significant discrepancy between the two analyses, it seems likely that the conversion we used was not quite correct.

The postfit residuals from the PEP and Tempo2 analyses (r_PEP and r_TEMPO2, respectively) are very similar. The corresponding two-sided Kolmogorov–Smirnov statistic is 0.035 (p-value 0.54), indicating that we cannot reject the null hypothesis that the two sets of residuals were drawn from the same distribution. We take as the nominal the average between these, r_A ≡ (r_PEP + r_TEMPO2)/2. A histogram of r_A is shown in Figure 1, where a Gaussian distribution is fit to the central part of the histogram yielding a standard deviation (SD) of 0.25 μs. By comparison, the SD of the full set of r_A is 0.57 μs. The figure shows that the "tails" of the data distribution are overpopulated, which is not unusual for real data. Figure 1 also shows a histogram of the difference of the two sets of residuals: r_D ≡ (r_PEP − r_TEMPO2). Because the residuals reported from Tempo2 are effectively reduced by the pulse offset, producing a set of values with an unweighted mean of zero, we have adjusted the residual difference by subtracting its weighted mean to make the two sets more directly comparable. The nearly rectangular distribution of r_D has a SD of 15 ns, which is 38 times smaller than the SD of r_A . The maximum absolute value among the r_D is 33 ns, making it unlikely that there exists a significant error in either of the two programs being compared. Finally, Figure 1 shows the lack of correlation between r_A and r_D and the calculated correlation coefficient ρ = 0.02 confirms that lack. The plot emphasizes that the residual differences r_D are small compared to the residuals r_PEP and r_TEMPO2, which is consistent with the two software packages having nearly identical, correctly programmed models.

Zoom In Zoom Out Reset image size

We have described the PEP and its principal capabilities. Over nearly six decades, this Fortran program of over 10⁵ lines has grown and evolved as scientific opportunities in the form of new data, usually of higher accuracy, and new data types have become available. That evolution has occasionally been facilitated by the direct comparison with independently written codes with which there is a substantial overlap in capability. We presented one such comparison (with Tempo2) here and mentioned another (with EPM). PEP is still in active use and the PEP source code is publicly available via GitHub and Zenodo. PEP is the only open source, general analysis tool for solar system astrometric data in the world, as far as we are aware.

We thank Ryan Shannon and George Hobbs for help with Tempo2. This material is based upon work supported by the National Science Foundation under grant No. PHY1708215, by the National Aeronautics and Space Administration under grant No. NNX15AC51G issued through the mission directorate, and under grant No. NNX16AH49H through the Massachusetts Space Grant Consortium.

The appendices describe the "Partial Pre-Reduction of the Normal Equations" (PPRNE), an unusual portion of the estimation module in PEP. It has the ability to transform (partially pre-reduce) the normal equations (NE) prior to obtaining a solution. The partially pre-reduced NE are smaller than the original but, for the parameters represented, yield the same parameter estimate and statistics as the full NE. Below, we start by discussing one motivation for having such a capability, Numerical Experiments. In Appendix B, Transformation of the NE, we provide a derivation of the algorithm and in Appendix C, Change of the Weighted Sum Squared Residuals, we show how to find the rms postfit weighted residual when PPRNE is used.

Several technical memoranda (TM) related to PEP are available in the public PEP documentation archive (Battat 2021a). This appendix is based on TM75-3 (19750729.pdf), which references TM74-1 (19740124.pdf). There are other TM available in the archive that describe algorithms in PEP that extend the work below. In particular, see TM75-2, TM80-1, TM88-06, and TM89-01.

A dynamic model of the solar system is limited in complexity by the ability of the available data to estimate the model-related parameters without making the NE degenerate (rank deficient) or even making the condition number unreasonably large. Thus, we are always using an "inadequate model" in our analysis, which implies that we must expect to get biased estimates. That is troubling in the case of those parameters that are of primary scientific interest. Many approaches have been proposed to address these concerns.

One approach involves running a large number of solutions that differ only in the inclusion or exclusion of members of a particular subset of the model parameters. The selection of solutions is a matter of judgment by the investigators.

The full ensemble of parameters included in a solar system solution often contains sets of "nuisance parameters" that must be included but whose estimates are of little or no interest. For example, when planetary radar observations were in regular use in the solar system analysis, we developed an ad hoc model of the topography of each planet in the region from which radar returns would come. These models (e.g., for Mars and Venus) could have hundreds of parameters each.

When the set of nuisance parameters becomes very large, the matrix inversion can take long enough that it interferes with the work, either because the cost or availability of machine time is limiting or because an investigator is trying to work in real time. As computers become more powerful and people tackle more difficult problems, this problem may become or cease to be important. We developed a solution dubbed "Partial Pre-Reduction of the Normal Equations." It entails a rigorous transformation of the NE to a new set of NE such that any solution of a subset of the new NE yields the same parameter estimate variance and covariance as would have been found for the parameters of the new NE by inverting the original NE. PPRNE is the central theme of these three appendices.

In classical WLS fitting, there is a data vector z of length m and a correct model H(X) where X is a vector of parameters of length n. No unique solution is possible for m < n. We define a loss function J = r^† R⁻¹ r, where r = z − H(X) is the prefit residual, R is the data noise covariance, and ^† indicates the transpose. When R is diagonal, the diagonal elements of R⁻¹ become a vector of data weights. The NE of WLS fitting take the form

$$\begin{eqnarray}&&B{\rm{\Delta }}=U,\mathrm{which}\ \mathrm{has}\ \mathrm{solution},{\rm{\Delta }}={B}^{-1}U,\end{eqnarray} \tag{ B1 }$$

where B is the n × n coefficient matrix, Δ is a vector of adjustments of length n, and U is the vector of length n that contains the error to be corrected. In particular:

$$\begin{eqnarray}\begin{array}{rcl}B & = & {A}^{\dagger }{R}^{-1}A\\ U & = & {A}^{\dagger }{R}^{-1}r\ \mathrm{where}\\ A & = & \displaystyle \frac{\partial H}{\partial X}.\end{array}\end{eqnarray} \tag{ B2 }$$

We divide the set of n parameters into two sets, α and β (of sizes n_α and n_β , respectively, with n_α + n_β = n), such that all of the "interesting parameters" are in the set α. This set includes all parameters whose estimate and uncertainties are of prime interest and all additional parameters that one sometimes may want to include in the analysis and other times may want to exclude from the analysis. The set β includes all nuisance parameters and all additional parameters that need to be part of the solution but can be ignored in the planned analysis.

We can order the parameters such that all α parameters precede all β parameters. Then the vectors and matrices of Equation (B1) can be written in partitioned form, where the upper left part is associated with the α set, the lower right part is associated with the β set, and the other parts have mixed association:

$$\begin{eqnarray}\begin{array}{rcl}B & = & \left(\begin{array}{cc}C & F\\ {F}^{\dagger } & D\end{array}\right),\quad {B}^{-1}=\left(\begin{array}{cc}P & Q\\ {Q}^{\dagger } & R\end{array}\right),\\ {\rm{\Delta }} & = & \left(\begin{array}{c}Y\\ Z\end{array}\right),\quad U=\left(\begin{array}{c}V\\ W\end{array}\right).\end{array}\end{eqnarray} \tag{ B3 }$$

From this it follows that

$$\begin{eqnarray}\begin{array}{rcl}V & = & {CY}+{FZ},\quad W={F}^{\dagger }Y+{DZ},\\ Z & = & {D}^{-1}(W-{F}^{\dagger }Y).\end{array}\end{eqnarray} \tag{ B4 }$$

By combining these we obtain

$$\begin{eqnarray}&&(C-{{FD}}^{-1}{F}^{\dagger })Y=V-{{FD}}^{-1}W,\end{eqnarray} \tag{ B5 }$$

$$\begin{eqnarray}&&Y={\bar{C}}^{-1}\bar{V},\end{eqnarray} \tag{ B6 }$$

where

$$\begin{eqnarray}&&\bar{C}=C-{{FD}}^{-1}{F}^{\dagger },\quad \bar{V}=V-{{FD}}^{-1}W.\end{eqnarray} \tag{ B7 }$$

Equation (B6) shows that the full NE can be replaced by a smaller set that yields the same adjustment for the parameters in the α set as would be obtained from the inversion of the full NE. We next show that the variance and covariance are also given correctly by the smaller NE

$$\begin{eqnarray}\begin{array}{rcl}{B}^{-1}B & = & \left(\begin{array}{cc}{PC}+{{QF}}^{\dagger } & {PF}+{QD}\\ {Q}^{\dagger }C+{{RF}}^{\dagger } & {Q}^{\dagger }F+{RD}\end{array}\right)\\ & = & \left(\begin{array}{cc}{I}_{\alpha } & 0\\ 0 & {I}_{\beta }\end{array}\right),\end{array}\end{eqnarray} \tag{ B8 }$$

where I_α and I_β are unit matrices of sizes n_α and n_β , respectively. To show that the PPRNE yields the same statistics as the original NE, we need to show that ${\bar{C}}^{-1}=P$. Then, using the result, PF+QD = 0, from Equation (B8):

$$\begin{eqnarray}\begin{array}{r}\begin{array}{rcl}P\bar{C} & = & P(C-{FD}^{-1}{F}^{\dagger })=PC-{PFD}^{-1}{F}^{\dagger }\\ & = & PC+{QDD}^{-1}{F}^{\dagger }=PC+{QF}^{\dagger }={I}_{\alpha }.\end{array}\end{array}\end{eqnarray} \tag{ B9 }$$

A metric in parameter space that yields a measure of the size N of the adjustment in the case of linear WLS fitting is the coefficient matrix B

$$\begin{eqnarray}&&{N}^{2}={{\rm{\Delta }}}^{\dagger }B{\rm{\Delta }}={{\rm{\Delta }}}^{\dagger }U.\end{eqnarray} \tag{ C1 }$$

For the predicted residual r_p , which is the postfit residual in the linear case,

$$\begin{eqnarray}\begin{array}{rcl}{r}_{p} & = & r-A{\rm{\Delta }}=r-A{\left({A}^{\dagger }{R}^{-1}A\right)}^{-1}{A}^{\dagger }{R}^{-1}r\\ & = & (I-G)r\end{array}\end{eqnarray} \tag{ C2 }$$

where $G=A{({A}^{\dagger }{R}^{-1}A)}^{-1}{A}^{\dagger }{R}^{-1}$ is a projection operator, G = GG. Then the sum squared weighted postfit residual is

$$\begin{eqnarray}&&{r}_{p}^{\dagger }{R}^{-1}{r}_{p}={r}^{\dagger }(I-{G}^{\dagger }){R}^{-1}(I-G)r.\end{eqnarray} \tag{ C3 }$$

We note that

$$\begin{eqnarray}&&{G}^{\dagger }{R}^{-1}={R}^{-1}A{({A}^{\dagger }{R}^{-1}A)}^{-1}{A}^{\dagger }{R}^{-1}={R}^{-1}G\end{eqnarray} \tag{ C4 }$$

and therefore

$$\begin{eqnarray}\begin{array}{rcl}{r}_{p}^{\dagger }{R}^{-1}{r}_{p} & = & {r}^{\dagger }{R}^{-1}(I-G)(I-G)r\\ & = & {r}^{\dagger }{R}^{-1}r-{r}^{\dagger }{R}^{-1}{Gr}.\end{array}\end{eqnarray} \tag{ C5 }$$

Revisiting Equation (C1) we find

$$\begin{eqnarray}&&{N}^{2}={{\rm{\Delta }}}^{\dagger }B{\rm{\Delta }}={r}^{\dagger }{R}^{-1}{Gr}\end{eqnarray} \tag{ C6 }$$

$$\begin{eqnarray}&&{r}^{\dagger }{R}^{-1}r-{r}_{p}^{\dagger }{R}^{-1}{r}_{p}={N}^{2}.\end{eqnarray} \tag{ C7 }$$

Thus, the effect of a fit is that the sum of the squared residuals drop (between prefit and postfit) by Δ^† U, which is easily calculated without first calculating the predicted residuals.

We next apply this analysis of residuals to the PPRNE. With the partitioning of the matrices and vectors into α and β parts, Equation (C1) becomes,

$$\begin{eqnarray}&&{N}^{2}={{\rm{\Delta }}}^{\dagger }U={Y}^{\dagger }\bar{V}+{N}_{\beta }^{2},\quad {N}_{\beta }^{2}={W}^{\dagger }{D}^{-1}W.\end{eqnarray} \tag{ C8 }$$

The quantity ${N}_{\beta }^{2}$ would be calculated once, say, at the time the NE were partially pre-reduced.