GipsyX/RTGx, A New Tool Set for Space Geodetic Operations and Research

GipsyX/RTGx is the Jet Propulsion Laboratory's (JPL) next generation software package for positioning, navigation, timing, and Earth science using measurements from three geodetic techniques: Global Navigation Satellite Systems (GNSS), Satellite Laser Ranging (SLR), and Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS); with Very Long Baseline Interferometry (VLBI) under development. The software facilitates combined estimation of geodetic and geophysical parameters using a Kalman filter approach on real or simulated data in both post-processing and in real-time. The estimated parameters include station coordinates and velocities, satellite orbits and clocks, Earth orientation, ionospheric and tropospheric delays. The software is also capable of full realization of a dynamic terrestrial reference through analysis and combination of time series of ground station coordinates. We present some key aspects of its new architecture, and describe some of its major applications, including Real-time orbit determination and ephemeris predictions in the U.S. Air Force Next Generation GPS Operational Control Segment (OCX), as well as in JPL's Global Differential GPS (GDGPS) System, supporting User Range Error (URE) of $<$ 5 cm RMS; precision post-processing GNSS orbit determination, including JPL's contributions to the International GNSS Service (IGS) with URE in the 2 cm RMS range; Precise point positioning (PPP) with ambiguity resolution, both statically and kinematically, for geodetic applications with 2 mm horizontal, and 6.5 mm vertical repeatability for static positioning; Operational orbit and clock determination for Low Earth Orbiting (LEO) satellites, such as NASA's Gravity Recovery and Climate Experiment (GRACE) mission with GRACE relative clock alignment at the 20 ps level.

Over the past thirty years, Global Navigation Satellite Systems (GNSS), starting with GPS, have become ubiquitous in our daily lives. A growing number of scientific, industrial, and public safety applications, spanning geodesy, geophysics, seismology, atmospheric sciences, weather forecast and climatology, time and frequency transfer, have come to depend on precision modeling of all aspects of GNSS at the cm and mm accuracy. Very few software packages exist with the capabilities to obtain these highest accuracy results, keeping up with the latest GNSS models and data processing techniques. Fewer still are capa-ble of also modeling and processing data from other space geodetic techniques such as SLR (Satellite laser ranging), DORIS (Doppler orbit determination and radio-positioning integrated by satellite), or VLBI (Very long baseline interferometry). Recently, JPL has embarked upon a major overhaul of its widely-used space geodetic data analysis software packages, GIPSY-OASIS and RTG, to derive GipsyX/RTGx, with many advanced new capabilities and quality features.
This paper is intended for the users and developers of space geodesy tools.
We present the design and some implementation notes for a new, but already widely used software set, and describe its performance in a broad range of applications spanning orbit determination and positioning operations, geodesy, and remote sensing. We first give some historical context, which influences any large software and research project design. We then present the overall design of the software, and sample use cases with tests of accuracy and precision.
The Jet Propulsion Laboratory (JPL) has a long history in space geodesy and precision orbit determination (POD) for Earth orbiting satellites, and Gip-syX/RTGx is the 4th major redesign of JPL's GNSS data analysis software.
For an overview of POD and measurement processing see the books [Tapley et al., 2004b] and [Bierman, 1977]. Born from Very Long Baseline Interferometry (VLBI) data analysis expertise and software [Sovers and Fanselow, 1987], the first JPL GPS data models were coded during 1984-1985 in the GPSOMC FOR-TRAN software [Sovers and Border, 1990]. Immediately thereafter we began development of a more comprehensive software set, the GPS Inferred Positioning System and Orbit Analysis Simulation Software package (GIPSY-OASIS, forever nicknamed GIPSY) [Wu et al., 1986], also coded in FORTRAN, to support the early campaign-style GPS orbit determination (GiG'91 for example [Melbourne et al., 1993]) and positioning experiments [Blewitt, 1989a]. That very early version of GIPSY should be considered a prototype for the more formally designed, enhanced, and documented version , written to support the first precision GPS flight experiment on the TOPEX/Poseidon altimetry mission [Bertiger et al., 1994] as well as the continuous GPS orbit determination operations at JPL that are the cornerstone for all precision GIPSY applications.
With the maturation of the GPS constellation, which became fully operational in 1994, real-time operational applications, such as differential positioning, started to emerge. The state-space formulation of real-time differential systems, today a universally accepted practice also known as wide area differential GPS (WADGPS) or global differential GPS (GDGPS), was conceived and patented at JPL [Yunck et al., 1989], and led to the creation of the Real Time GIPSY (RTG) software. RTG ported most of GIPSY's precise satellite and signal models from FORTRAN into C, but redesigned the software architecture, and in particular the estimation filter, for efficient real-time operations.
Successful early demonstration of RTG [Whitehead et al., 1998] has led the United States Federal Aviation Authority's (FAA) to adopt JPL's state-space algorithms for its Wide Area Augmentation System (WAAS), with RTG as the prototype software for this critical aviation infrastructure, which has since been adopted by several other countries as well [Bertiger et al., 1997]. In 2000, RTG became core software for JPL's GDGPS system, as well as for other commercial systems.
As with any JPL-developed software, GIPSY and RTG are owned by the Cal- project to develop and deploy OCX [http://insidegnss.com/raytheon-wins-1-5-billion-gps-ocx-contract/], [Bertiger et al., 2010a], and we set out to create JPL's 4th generation GPS data processing software, entitled RTGx. Later, the NASA Space Geodesy Project (SGP) provided additional support to enhance post-processing for space geodesy [Merkowitz et al., 2018]. All programs synergistically benefited.
This was an opportunity to modernize the entire geodetic data processing software set at JPL and unify the agile real-time processing of RTG with the GipsyX/RTGx, as was GIPSY-OASIS, is a member of a small 'club' of available precision space geodetic data analysis packages, with diverse capabilities, and with a track record of contributing high quality GNSS data analysis products. These include Bernese [Dach et al., 2015], EPOS.P8 [Uhlemann et al., 2015], GAMIT [Herring and Floyd, 2018], GEODYN [Pavlis et al., 2017], GINS [Marty et al., 2011], NAPEOS , and PANDA [Shi et al., 2010].  • Hot start capability with archived filter state and covariance • Integer GNSS phase ambiguity resolution in post-processing and in realtime • Inertial measurement modeling to optimally combine GNSS and inertial data in real-time and in post-processing • Hyperlinked documentation using Doxygen and Sphinx • User input interface featuring inheritance • Reference frame determination and time series analysis of Earth station positions Table 2 shows the size of the entire GipsyX/RTGx distribution in terms of the number of lines of code. Raw lines are just the total number of lines in all the source files. Physical lines are lines that are not comments or blank. The majority of the C++ code occurs in the single executable module, "rtgx", which models the GNSS data and other tracking data types, and estimates parameters such as station and satellite position. We describe the structure of "rtgx" below.
There are approximately 140 individual executable programs in GipsyX/RTGx.
Included are many small utilities, such as time format manipulation, time series and data manipulation, trend analysis, data archive access and plotting utilities.
Most users will only use a small set of these. In addition to the tracking data analysis tools, there is a set of Python3 executables dedicated to time series analysis of positions on the Earth and reference frame determination as well as many other small utility tools (for example, to convert between file formats or timing formats).
We have attempted to apply in this development the many lessons learned from the previous iterations of GIPSY and RTG. GIPSY was written by many people under a fairly loose set of coding rules. The main computational pieces of GIPSY were divided into several independent FORTRAN modules (orbit integration, signal/Earth model, filter, smoother...). There was considerable inconsistency in the input to these pieces, which were mostly FORTRAN namelists.  [Bertiger et al., 2010c]. GipsyX/RTGx improves upon GIPSY in two ways here. One, there is a uniform input format, referred to as a 'tree' input defined in detail below, which has a general structure similar to YAML. We did not use YAML because, at the time of our development, it did not have some object-oriented inheritance properties or arbitrary script execution that would allow human readable, compact input for GNSS constellations where many of the satellites have identical properties. Second, instead of writing many individual executable modules as in GIPSY, the object-oriented C++ allowed us to write a single main executable, rtgx.
For version control, we are using svn, subversion [Pilato et al., 2008] chosen before distributed version control systems such as git [Chacon and Straub, 2014] were more popular. We believe that documentation is best maintained inside the source code; thus, we chose to use Doxygen for C++ code, and Sphinx for Python. Doxygen and Sphinx both output html documentation and it is easy to link the two sets of documentation. The mathematical description for the radiometric signals, solid tides, Earth orbiting satellite forces, ordinary differential equation integrator, etc., are best written in L A T E X. Our html documentation includes a pdf containing the mathematical description, generated with L A T E X, with links to the source code documentation.
Having a single main executable, nearly all use cases can be constructed as a very thin wrapper around a single command line: rtgx myInput.tree To automate the most common tasks carried out by our user community, namely precise static point positioning with single receiver ambiguity resolution, we have provided a Python3 wrapper, where the only required input is a RINEX2 (Receiver Independent Exchange Format) or RINEX3 GNSS data file.

gd2e.py -rnx my.rnx
This is the analog of GIPSY's gd2p.pl automation script, and it contains many enhancements over its predecessor in terms of ease of use and documentation.

User interface -the input tree
Several of our executables, especially those requiring more complicated and finely controllable user inputs, are configured by a single interface called a tree.
A tree is a text file with a hierarchical structure identified using Python-like indentations that consist of roots (the highest level), branches, and leaves (the lowest level; the leaves of the tree are the parts of a tree with no branches and are generally the ones that contain specific data). See electronic supplement 1 for details.

Main C++ Software Modules/Classes
The computationally intensive software breaks into several broad categories: • Data editing • Orbit Integration and force models All of these items are implemented as C++ classes or functions in a single executable module, rtgx, except for data editing, which is provided by the "gde" module.

gde Data Editor
GNSS data must be edited for phase breaks and gross outliers before it is processed to fit model parameters. The GNSS Data Editor (gde), is a C++ code that implements two distinct algorithms: turbo-edit [Blewitt, 1990], and simple continuity checks of linear combinations of phase data based on lowdegree polynominal fits. Turbo-edit looks at averages of pseudorange − phase measurements to detect jumps in the phase measurements. It is well suited to low rate data, such as the large number of RINEX files typically recorded at a 30-sec rate. When data rates are high, for example sampled every second as often seen in real-time applications, it is possible to detect phase breaks simply by monitoring the change in time of differences in phase observations on two frequencies. Such a difference cancels out receiver and transmitter clock jumps.
Removing a low order time polynomial from these differences effectively removes the slow-changing ionospheric signal in this difference. Single differences of dualfrequency combinations across two satellites are also available as an option.

Orbit Integration and force models
To model the path of a satellite in orbit about the Earth, we implement precise force models that include gravitational forces and surface forces on the satellite. These forces may contain model parameters. Given initial values for the satellite epoch state (position and velocity at the initial epoch time) and parameter values, the second order differential system of equations for the time evolution of the satellite state and the partial derivatives of the satellite epoch state with respect to the model parameters are integrated numerically forward in time. Force models represented in GipsyX/RTGx include the Earth's gravity field as a spherical harmonic expansion to arbitrary degree and order including time variability of these coefficients, gravitational effects of Earth solid tides, pole tides and ocean tides, first order General Relativistic acceleration and atmospheric drag (DTM-2000, [Bruinsma et al., 2003]). Conical models for solar eclipse by the Earth and Moon are implemented. Generic Earth albedo ( [Knocke, 1989;Knocke et al., 1988]) and solar radiation pressure ( ) models are available through a flexible tree input interface which enables the specification of arbitrary sets of panels along with their orientation and reflective properties. More specific satellite radiation pressure models, including some empirical GPS solar radiation pressure models specially developed at JPL ([Bar-Sever and Kuang, 2005] and [Sibthorpe et al., 2011]), may also be input either as arbitrary Fourier series, or bi-linearly interpolated tables dependent on Sun azimuth and elevation angles in the satellite body-fixed frame.
In order to accurately compute surface force models, spacecraft orientation is required, and may be input in a variety of ways, including time-series of attitude quaternions and, especially useful for simulations, dynamic pointing specified via mathematical operations on vectors from the satellite to the Sun and Earth, and the satellite's velocity with a defined input calculus. Attitude may also be selected from a number of specific models including those for GPS ([Bar-Sever, 1996;Bar-Sever et al., 1996]), GLONASS ( ), Galileo (https://www.gsc-europa.eu/support-to-developers/galileo-iovsatellite-metadata) and configurable yaw steering and orbit normal modes for BeiDou geostationary and inclined geostationary satellites ( [Hugentobler and Montenbruck, 2017]). For other orbiting satellites including altimetry missions, attitude may be selected from models or quaternions.
We have implemented a general-purpose ordinary differential equation (ODE) solver which we use to obtain numerical solutions of the time evolution of the satellite state. This ODE solver includes multistep methods [Hairer et al., 1993] as well as high-order embedded Runge-Kutta methods [Prince and Dormand, 1981]. The solver employs typical local-truncation-error-estimate adaptive step size techniques. In order to allow for discontinuities in the integrated source functions (e.g., the solar radiation pressure transition between shadow regimes, and satellite thruster firings have discontinuous time derivatives) and yet still employ high-order integration techniques, the ODE solver employs a modelspecific discontinuity detection logic, following the"gstop-function" approach specialized at JPL by Fred Krogh for integrating spacecraft trajectories ( [Krogh, 1972]) https://mathalacarte.com/fkrogh/pub/dxrkmeth.pdf.

Signal Model
The basic signal model is the time of propagation from a transmitter to a receiver. For GNSS range measurements recorded at receiver time,t r , determined by the receiver's clock and transmitted at timet t , we model the measured range where c is the speed of light, d trop is a delay due to troposphere and d iono is the delay due to the ionosphere. For many GNSS measurements the first order delay of the ionosphere is removed by a dual-frequency combination of the data ] and a model for higher order effects [Kedar et al., 2003].
The troposphere is modeled as a delay at zenith plus gradient parameters ([ Bar-Sever et al., 1998], [Böhm and Schuh, 2004], [Böhm et al., 2006], [Böhm et al., 2015]). The difference in the time of reception,t r , and the time of transmission, t t , in equation 1 can be modeled as wheret r −t r is the difference between the time on the receiver clock and proper time(time clock would have read with no errors) andt r − t r is the difference between proper time and coordinate time (General Relativistic effects, [Moyer, 1971], [Thomas, 1975], [Moyer, 1981]). The last two lines of eq. 2 contain similar differences for the transmitter. The middle term, difference in coordinate time, For multi-GNSS, eq. 1, must be modified due to different delays of ranging codes through the receiver for different constellations. The code adds a bias parameter for each receiver by constellation [Odijk and Teunissen, 2013]. One constellation must be held as reference to prevent singularity. Typically we hold GPS as reference, but it is arbitrary.
The model for GNSS phase data is identical to range, except two terms must be added to eq. 1 where φ is the modeled phase, B t r , is an arbitrary phase bias between the transmitter and receiver over time periods where the receiver has not lost lock, and ω(t r ) is the phase windup [Wu et al., 1993]. The B t r may further be modeled as integer number of wavelengths and fractional hardware delays in the given receiver and transmitter; see [Blewitt, 1989b] for treatment of network ambiguity resolution and [Bertiger et al., 2010c] for treatment of the integer ambiguities with a single receiver.

Earth Models
For receivers or transmitters that are not orbiting the Earth and are either attached to or associated with the Earth's crust, we must model the deformation of the Earth's crust and the orientation of the Earth in inertial space since satellite positions are integrated in inertial space. GipsyX/RTGx implements the IERS standards [Petit and Luzum, 2010]. Adjustable parameters include the Earth orientation parameters, polar motion and hour angle.

Filter, Smoother, Ambiguity Resolution
Except for single receiver use cases, most of the computation time is dominated by filtering, smoothing and ambiguity resolution functions. If we treated all the parameters as constant in time, the filter and smoother parts of the code would just be a standard least squares; but many of the parameters do vary in time, for instance the zenith troposphere delay or the receiver and transmitter clock errors. Many force parameters, including empirical accelerations affecting satellite dynamics, are also best treated as stochastic variables in a process called reduced-dynamic orbit determination [Yunck et al., 1990;Wu et al., 1991].
We chose to handle these time variations as a first order Markov process with a forward-in-time filter implemented as a Square Root Information Filter (SRIF) detailed in Bierman [Bierman, 1977]. The square root implementation allows for better conditioned matrices, and the use of Householder (orthogonal transformation in multi-dimensional space) transformations in our implementation of the Bierman algorithms lends itself naturally to more efficient use of processor specific pipelining [Jeong et al., 2012], when available, along with the Message Passing Interface (MPI [Gropp et al., 1999]) for processing across multiple CPU cores and computer clusters (Beowulf [Sterling et al., 2002]).
The forward filter is the best fit of all the past data up to the current epoch.
In order to have the parameters fit the future as well as the past data optimally, one needs to smooth the information back in time. We again follow Bierman's algorithm [Bierman, 1977] for a SRIF smoother, which uses a series of Givens (two-dimensional orthogonal rotation) and Householder transformations. We again optimize efficiency using machine-specific-low-level routines and MPI where appropriate. After editing for outliers in the smoother, integer ambiguities may be constrained for GNSS data. For single receiver ambiguity resolution, we follow the procedure in [Bertiger et al., 2010c] extended to multi-GNSS, applying constraints to single differences. For multiple GNSS receivers, we can resolve double difference integer ambiguities following [Blewitt, 1989a] extended to multi-GNSS to constrain the double-differenced phase ambiguities, as well as adding single-receiver constraints when appropriate external information is available. The single receiver methods are related but not identical to methods developed in [Laurichesse et al., 2009;Laurichesse, 2011], [Loyer et al., 2012]. Our ambiguity resolution algorithms, do not currently extend to GLONASS Frequency Division Multiple Access (FDMA) signals.

Main executable, rtgx
A detailed description of the main loop of the primary executable is contained in electronic supplement 2, rtgx main loop.  [Johnston et al., 2017], two sets of GPS products characterized by different latencies, precision and accuracy levels, but all based on filtering and smoothing of ground tracking data in post-processing."Rapid" products are generated and distributed daily. "Final" products are generated weekly and typically have a 14-day latency. We also distribute "Ultra-rapid" products, generated every hour. While ultra-rapid and rapid products are tightly tied to the IGS realization of the ITRF by fixing the coordinates of a subnetwork of stations in the POD process to their ITRF coordinates, Final products come in three flavors. Fiducial-free products, identified by the suffix 'nf' (non-fiducial) in their names, are not tied to any specific terrestrial frame; the positions of the stations constituting the network are left floating with a 1 km a priori sigma, so that the frame of these solutions changes from day-to-day.
No-net-rotation ('nnr' suffix) products are such that 3 no-net-rotation (relative to the ITRF) constraints are enforced in their generation process. Final products with no suffix are tied to the ITRF by means of 7 constraints (3 rotations, 3 translations and 1 scale parameter) and are referred to as NNRTS (no-netrotation, translation, or scale) or 'fiducial' products, though their fiducialization differs from Rapid and Ultra-Rapid. Every set of products contain files with consistent orbital state estimates, transmitter clock estimates, spacecraft attitude information, Earth rotation parameter (ERP) estimates, and widelane phase biases information. In addition, the no-net-rotation and fiducial-free solutions are associated with coordinate transformation files referred to as x-files. These x-files contain the set of 7 Helmert parameters (3 small rotations, 3 translations, and 1 scale parameter) needed to transform from the frame of the day to the ITRF. The format of these different products is detailed in the GipsyX software documentation. Following extensive testing and validation, the JPL GNSS analysis center transitioned seamlessly from products generated using the GIPSY software package to products created using GipsyX on January 29, 2017. More recently, the entire span of 1994 through 2018 were reprocessed using GipsyX, consistent with IGS repro2 standards. All products are available at https://sideshow.jpl.nasa.gov/pub/JPL_GNSS_Products. Table 3 displays the characteristics of each type of post-processed product.
Latencies are indicated as well as commonly used performance metrics. JPL orbit and clock products cover 30 hours centered at noon, so that each daily solution overlaps with the next-day solution by 6 hours. After removing 30minute tails at both ends of this overlap period, RMS statistics on orbit and clock differences are computed in the 5-hour overlap period as an internal measure of the precision of the product. The long-term medians of the RMS clock and orbit differences in the overlap is shown as "Precision" in Table 3. We measure accuracy in Table 3 by comparing JPL's orbit and clock solutions with the combined final IGS orbit and clock solutions for the same days. Although we have listed this as "Accuracy" in Table 3, we note that the IGS Final product is a weighted average of all the contributing analysis centers including JPL.    . Orbit modeling differences, e.g., in solar radiation pressure modeling, are known to impact the accuracy of GNSS-based ERP determination. In particular, they could explain the bias visible in Fig. 3 [Bizouard et al., 2017]. Table 4 of [Bizouard et al., 2017] shows Xp and Yp scatter of IGS combined at 31 and 27 µas respectively and 10 µs for LOD for ITRF 2014 from 2010-2015.

Four Constellation GNSS Orbits and Clocks
A greater number of satellites may provide significant benefits to users, including additional coverage in urban canyons, and perhaps better estimates of parameters such as troposphere, Earth orientation, and geocenter. To this end, many IGS Multi-GNSS Experiment (MGEX, , http://mgex.igs.org/) analysis centers are now routinely processing multiple   to 2.6 mm at the surface of the Earth over a day, is attributed to orbit modeling errors and in particular solar radiation pressure modeling errors in GNSS data processing, as discussed in [Sibthorpe et al., 2011] and [Ray, 1996] for instance.
constellations, particularly GPS, BeiDou, Galileo and GLONASS, although as yet no official combined product is available. It is worth considering however, that processing new satellites might also have unintended impacts on, e.g. GPS quality compared to existing GPS-only products, because such new satellites may require individual treatment to allow for accurate modeling of their specific observation types (e.g. GLONASS frequency dependent code biases) and the spacecraft themselves (e.g. solar radiation pressure, attitude and antenna offsets [Montenbruck et al., 2015]). Although JPL has been processing these constellations for some years now in real-time (see section 6.2), we are still refining our station-data selection procedures for fully automated post-processed products, as well as our force model choices, alongside internal quality testing, before making an operational version corresponding to one or more of our GPS-only product lines (see section 6.1.1). Such internal testing showed a dra-

Real-time GNSS Orbit and Clock solutions, GNSS Differential Corrections
The real-time configuration of GipsyX/RTGx is typically called just RTGx.
RTGx has been driving real-time orbit determination processes in the GDGPS System since 2014 (RTG was the engine prior to that, with some overlap). These For the full derivation, see the appendix in [Zumberge and Bertiger, 1996].
Here we have used the standard approximations of the coefficients. The realtime URE relative to post-processed solutions is typically 5 cm RMS (Fig.   4).
The specialized requirements of real-time operations have driven the design of some unique and powerful RTGx features. One such feature is the ability to reconfigure a real-time running filter into partitions such that some parameters associated with specific satellites or specific stations, are estimated but have no influence on any other parameters. At the GDGPS System these "decoupled" partitions are used to accommodate unhealthy satellites in order to keep monitoring them while protecting the rest of the estimated constellations from the potential mismodeling of unhealthy satellites, for example due to maneuvers.
Our decoupling algorithm is equivalent to infinitely de-weighting the data from the associated satellite or station following a white-noise re-set of the parameters associated with the satellite or station. For earthquake monitoring the GDGPS System operates orbit determination filters with multiple decoupled partitions where the position of hundreds of ground sites (a site per partition) are estimated kinematically at 1 Hz [Larson et al., 2007] , but these sites are  Table 3 above for post-processing accuracy). The RMS URE approaches 5 cm with ambiguity resolution and 90 or more real-time tracking sites. A cold start of a real-time GNSS orbit determination may take from hours to days to fully converge (i.e. to achieve the accuracy level of a long running continuous filter), depending on the quality of the initial states. In certain operational scenarios it is desirable to restart the filter with some modified configuration, while achieving instant convergence. This capability is enabled in RTGx through the commanded (or scheduled) saving of a "snapshot" of the entire filter state and model state, and the ability to read this snapshot upon restarting to achieve a hot-start capability. By instant convergence we mean that if the identical data are fed to the filter after the snapshot epoch, the solution will be identical to the continuous running filter. This feature allows the operator to rapidly correct for errors in the past and catch back up to real-time.
Other specialized features designed to support robust real-time operations include special handling for clock references to ensure smooth transition from one reference clock to another, on-the-fly decoupling of satellites or sites, pairwise filter operations combining a low cadence orbit filter (e.g., every 60 seconds) with a high cadence (e.g. every second) clock filter for optimal throughput.
The GDGPS System produces a variety of real-time differential corrections for the GNSS broadcast ephemerides, based on the real-time orbit and clock solutions, including formats that enable real-time ambiguity resolutions such as the state-space representation standards 10403.3 of the Radio Technical   , and are selected to be well-distributed geographically. Figure 6 shows the geographic distribution of sites used for this test. First, each site was point positioned using a tree very similar to the default tree included with GipsyX-1.0. The only substantive difference with the default tree is the use of Ionosphere Exchange (IONEX) files [Schaer et al., 1998] [Böhm et al., 2015] nominal troposphere and mapping function was used.
Each site was positioned using JPL's repro3.0 IGS14 Final products, with separate PPP runs for the non-fiducial (NF), no-net rotation (NNR), and no-net translation, rotation, or scale constraint (NNRTS also referred to as "fiducial") products from 2008-06-01 through 2008-11-30. In the case of the NNR and NF products, a 7-parameter Helmert transformation is performed using the appropriate product x-file to convert the resulting position into the reference frame.
Positions are compared with those given in ITRF2014. For each station, a 5sigma edit is performed on total position differences, then the root-mean-square (RMS) is calculated in east, north and vertical (E, N, V) components. We then report the median of each of these components in Table 5. Figure 7 shows a histogram of the frame repeatability for all 59 stations including all outliers.
All three methods (NF, NNR, NNRTS) produce nearly identical results. When compared to the GIPSY-OASIS results in [Bertiger et al., 2010c], E and N are almost identical but the results here are a bit worse in V (6.5 vs. 6.0 mm).
We attribute the difference in vertical to station selection differences (i.e. 59 IGS14 stations vs. 106 IGS08), as vertical repeatability is highly site-dependent.
While using IONEX files does not improve repeatability in this analysis, we still recommend its use to be consistent with the models used in the Final POD process, and because other analyses show that omitting it can introduce an offset in the z-component [Ries et al., 2018]. In this test, the offset is an average of 0.4 mm. The ENV frame repeatability of about 2, 2, and 6 mm in Table  5 is significantly better than the 24-hr RMS repeatability of 3, 6, and 9 mm shown in [Soycan, 2011] using Bernese 5.0 software. We note that significant improvements are seen in the east component with bias fixing [Bertiger et al., 2010c] with the GipsyX/RTGx bias fixing methods as well as the result with PANDA software [Ge et al., 2007] with bias fixing based on fractional satellite delays. In [Ge et al., 2007], they compute repeatability relative to IGS weekly solutions after removing a 7-parameter Helmert transformation and get mean RMS E, N, and V repeatability of 2.4, 2.7, and 5.3 mm with bias fixing. These are not exactly the same statistics as shown in Table 5. The weekly solutions will take out some vertical signal from seasonal loading not contained in the ITRF2014 frame. Figure 7: Histogram: East, North, and Vertical ITRF repeatability for 59 sites used for static PPP testing using NF products and IONEX corrections. Each count represents a single day of processing for a single station. Figure 8 shows the seasonality of PPP residuals relative to the ITRF solution for ALGO (Algonquin Park, Canada) using IONEX corrections and NF products for 16 years years. Most of the vertical residual can be explained by atmospheric and hydrological loading [Tregoning and Watson, 2009], models of which are also plotted. The horizontal displacements also show seasonality due to other effects such as snow cover [Larson and Nievinski, 2013]. This example illustrates one complexity of repeatability calculations, and how results can be biased by short time periods and/or inhomogenous distributions of stations.

Kinematic Positioning
For many applications the receiver is moving, and it is the accuracy of the kinematic position of the GNSS receiver that matters. Applications, among many, include synthetic aperture radar (SAR) from aircraft, ocean floating buoys used to determine sea surface height, co-seismic deformation, and ice sheet movement, see [Hensley et al., 2008;Bonnefond et al., 2003;Born et al., 1994;Simons et al., 2011;Doake et al., 2002] for example. To give an idea of the expected performance from GipsyX, we kinematically positioned a set of 45 ITRF2014 frame stations every 5 minutes for the month of June 2008 using GPS data. Table 6 shows the RMS differences with the ITRF2014 frame positions where the position has been adjusted as a loose random walk, 1m/ √ s, while the troposphere is a bit more constrained than in the static point positioning. Instead of adjusting the zenith delay and tropospheric gradient parameters, only a more constrained zenith delay is adjusted to help break the correlation with the random walk vertical position. and the dry zenith delay may be computed using a pressure sensor. To demonstrate the expected performance of kinematic positioning with better known troposphere, we repeated the experiment from Table 6, but fixed the zenith delay and the gradient parameters to the values from static 24-hour position solutions. Comparing Table 6 with Table 7, one can see a significant and larger improvement in the vertical.

Precise orbit and clock determination of low Earth orbiters
Both GIPSY and RTG supported POD of satellites in Low Earth Orbit (LEO), starting with the launch of TOPEX/Poseidon [Bertiger et al., 1994].
As a superset of these two software, GipsyX/RTGx supports precise postprocessing, reduced-dynamic orbit determination, as well as embedded real-time POD. For post-processed LEO POD using GPS data, the GipsyX/RTGx provides identical capability to GIPSY, but with a simpler user interface along with additional capability for the newer GNSS constellations, DORIS, and SLR. Tests were performed with GRACE [Tapley et al., 2004a], the followon GRACE mission (GRACE-FO), and Jason-2 [Lambin et al., 2010;Bertiger et al., 2010b;Cerri et al., 2010]. Since the GPS determined orbits for Jason-2 are sub-centimeter accurate [Bertiger et al., 2010b], we can use those orbits to measure our accuracy for Jason-2 orbits determined with other data types.
Mixing GNSS, SLR, and DORIS data is also possible in GipsyX/RTGx, but has not had extensive testing yet. In the subsections below, we detail the expected accuracy and software capabilities for LEO POD using GPS, SLR, and DORIS data.

SLR Orbit determination
GipsyX has many capabilities for SLR processing, including producing residuals to orbits determined from other techniques and independent orbit determination. For one experiment, We used 30-hour arcs of SLR data from 15 well-distributed, low-bias SLR sites, centered on noon of each day in 2015 to de-termine Jason 2's orbit with a median radial RMS difference to GPS-determined orbits of 1.9 cm. We also validated our SLR capabilities by performing orbit determination on the LAGEOS 1 and 2 satellites using 7-day arcs for all of 2017. We then compared our orbits to the ILRSA combined orbits. Our LA-GEOS 1 and 2 orbits have an overall RMS difference of 5 mm in radial, 16 mm in cross-track and 18 mm in along-track to LAGEOS, which is comparable in magnitude to the inter-center differences measured by the ILRSA combination to other centers over the course of the year (averages of 5 mm, 24 mm, and 25 mm in the weekly ilrsa weekly sum files across all centers for both LAGEOS satellites).

GPS -GRACE, GRACE-FO POD and time synchronization
GipsyX/RTGx is the operational software used for GRACE-FO to determine the orbit and clock of the GRACE-FO spacecraft.
GRACE and GRACE-FO are among the best low Earth orbiting satellites to test POD with GPS since they have a measurement of the inter-satellite range (up to a bias) that is good to the micron-level. This is the dual-one way range measurement made at K and Ka band, K-Band Range (KBR) [Dunn et al., 2003;Thomas, 1999]. It is the KBR's precision along with an accelerometer measuring the non-gravitational forces that allows the precise recovery of the Earth's time varying gravity field. Although GipsyX/RTGx, can process KBR range data and use accelerometer data, neither is used in the operational POD for GRACE or GRACE-FO discussed here. Figure 9 shows the daily standard deviation of the difference of the range determined by GPS POD and the KBR (clocks aligned with the independent operational code). The processing used GPS data sampled every 10-sec, antenna calibrations based on one month of data, August 2010, reduced dynamics selected to perform well over the full time period, and bias fixing. The average daily baseline standard deviation was 1.6 mm. This is significantly better than the 10 mm found in Table 8 of [Kang et al., 2006] using the Center for Space Research's MSODP software which does not include bias fixing. It improves over the results in [Bertiger et al., 2010c] mostly due to the use of higher rate data. [Kroes et al., 2005] and [Jggi et al., 2007] show improved baseline determination of a little under a mm if biases are fixed between the two spacecraft with a simultaneous POD of both spacecraft.
In this processing, both spacecraft are processed independently. Although the KBR measurement is not strongly dependent on the GPS POD/Clock solution for GRACE at the 100 micron level, it is necessary to have the error of the relative clocks between the two GRACE spacecraft less than 160 ps, 1-sigma, after removing a bias (see [Thomas, 1999] eq. 3.19) to have the errors in the KBR range due to the clock alignment less than 0.5 microns. Figure 10 shows a measure of the precision of the relative clock between the two GRACE spacecraft. Since we are solving for the GRACE clock as a white noise process, relative to a fixed GPS reference with 30-hours of data centered on noon of each day, there is a six hour overlap from one day to the next. If our reference clock were the same, and our solution were perfect, the difference of the clock solutions between these two days would be zero. In eq. 5, the first term is the difference of the two GRACEA clocks, the second GRACEB. Taking the difference of the A and B differences removes any reference clock, thus the RMS of this quantity over the central 5 hours of the 6-hours is a good measure of the relative clock precision. Fig. 10 shows that we are significantly under the 160 ps requirement with a median value of 20.3 ps. In addition to GNSS phase and range data, GipsyX/RTGx can process almost any radio-metric phase and range data and is easily extended to data types that are linear combinations of radio-metric range/phase data; for instance, Satellite Laser Ranging or Doppler Orbitography and Radiopositioning Integrated by Satellite (DORIS) . Current DORIS data is given in a RINEX like format [Auriol and Tourain, 2010] with dual-frequency phase measurements and range measurements similar to GNSS. Traditionally this is processed as Doppler where φ(τ ) is the measured phase on an orbiter from a ground transmitter at it's local time τ . Prior to June 20, 2008 the Doppler measurements were given directly with corrections for errors in the satellite clock and phase center offsets.
The count time, τ 1 − τ 2 , is 10 seconds of local satellite time in eq. 6. When forming DORIS Doppler from the phase on the RINEX DORIS data, one must model the satellite receiver clock errors for the computed Doppler model, C, The phase measurement, with units of length, is written as φ(τ ) = ψ(τ )+cE(τ ), where E(τ ) is the difference between coordinate time (t) and the local time of the receiver in seconds and c is the speed of light. For the DORIS system, the receiver is driven by a stable clock (Allan deviation on the order of 10 −13 at a day) and the rate of the clock errors including relativistic effects, E(τ1)−E(τ2) t1−t2 is small, typically on the order of 10 −7 s/s, we can expand eq. 7 The computed measurement model in eq. 8 is very close to the difference of phase measurements and fairly easy to code in terms of the basic phase measurement used for GNSS. This model also may be used for the legacy DORIS doppler data since the clock error terms are essentially zero, due to pre-calibration. On the DORIS RINEX file, the nominal receiver clock values are given as a time series derived from the range data. A quadratic fit over time periods on the order of a day maybe used for E(τ ).
We have tested this with data from the Jason-2 spacecraft [Couhert et al., 2015] using station coordinates from the latest DPOD (DORIS Precise Orbit Determination) solution, DPOD2014, aligned on ITRF2014 [Moreaux et al., 2019]. To take into account the effect of the South Atlantic Anomaly (SAA) on the onboard oscillator [Willis et al., 2016], we did not use any correction model but removed a few stations in the South America region. In Fig. 11, we compare the Jason-2 radial orbit position determined with the DORIS RINEX data from Feb. 14, 2014 through August 23, 2014 with a GPS-determined orbit whose radial accuracy is about 5 mm [Desai et al., 2018] using ITRF2014 coordinates for tracking stations and typical RMS differences with the GipsyX/RTGx DORIS determined orbit are at the 9.9 mm level. Errors in radial position go directly into the mission's prime objective to measure sea surface height and are the most important metric for ocean altimetry missions.
To better evaluate the SAA effect on Jason2, we will soon investigate which stations are the most affected and use the new available correction model from [Belli and Exertier, 2018], based on T2L2 measurements. In the future, we plan to analyze the Sentinel-3A DORIS RINEX data. The fact that satellite clock for Sentinel-3A is the same for the GPS and the DORIS on-board instruments should allow us to compare the GPS clock solution and the DORIS clock. The GPS solution for the on-board clock could be used to properly correct the SAA effects for this satellite, [Jalabert and Mercier, 2018].

Multi-Technique Geodesy
As an example of multi-technique use, we processed GPS and SLR data simultaneously using GPS and SLR data from low Earth orbiters to tie the systems together. For the development of reference frames, the ties are traditionally made through local ground surveys to measure the vector between the GPS and SLR instruments. With only three days of data, and no ground survey the GipsyX/RTGx solution reproduced these independent ground ties to 1 cm (3D). The tie between SLR and GPS in the GipsyX/RTGx solution comes solely from the space ties on Jason-1 and the GRACE tandem, testifying not only to the benefits of combining data at the observation level but also to the fidelity of the GipsyX/RTGx modeling.

Building Terrestrial Reference Frames and Fitting Ground Site Time Series
GipsyX includes station coordinate processing software which can combine solutions from several techniques to build reference frames, simulate future performance, compute transformation parameters between reference frames, and fit individual time series to estimate positions, velocities, seasonals, and possible discontinuities in the linear and periodic fit, which we refer to as breaks.
Each tool has command line help available. On-line training provides detailed instructions and example data sets for the most common use cases.    in Table 8 and DTRF2014 [Seitz et al., 2016] in Table 9.
See electronic supplement 3 for details.
GipsyX includes a network simulation tool which can generate output files at a given frequency for a specified time span, for example daily files spanning one year. The user provides a model consisting of positions, velocities, breaks, and/or seasonal terms. Various forms of white noise can then be injected. Coordinate noise can be added which is independent from site to site. Geocenter, scale, and/or rotational noise can be injected which is correlated and impacts all sites. Simulated solutions can be processed just like real ones to test software or predict performance of future networks.
The transformation between two reference frames at a particular time can be described by three translations, three rotations, and one scale. Daily transformation parameters known as x-files are computed using an input reference frame to predict coordinates at a particular epoch and then using those predictions to estimate transformation parameters to GNSS positions observed at that epoch. Daily x-files are publicly available as described in section 6.1.1. A plot of x-file parameters between IGS14 and GipsyX free-network solutions is shown in Figure 13 . The TX parameter has a mean value of -0.5 mm and a standard deviation of 4.7 mm. The TY parameter has a mean value of -0.7 mm and a standard deviation of 5.6 mm. The TZ parameter has a mean value of  shown that adding LEO satellites such as GRACE and Jason allow GNSS to be competitive with SLR for both linear and annual geocenter estimation .

Summary
The new GipsyX/RTGx software is a robust and powerful tool for geodetic data analysis and simulations, incorporating decades of expertise and lessons learned from the design and operations of previous software generations, and their applications to the most challenging positioning, navigation, timing, and science applications. GipsyX/RTGx now underlie all GNSS orbit determination operations at JPL, and with hundreds of academic and research licenses, powers geodetic analyses and science operations across the globe.