In general, GPSdancer strictly follows the IERS 2010 conventions and the processing standards of the IGS REPRO 2 campaign.

As GPSdancer was written entirely in JAVA, anybody who wants to have (...tested!) JAVA implementations of the models below can contact the project. Bizarly, the world of science is still full of people who refuse to consider anything else than FORTRAN, even claiming that it is somehow very efficient. Please note that the FORTRAN compilers that you use today were developed for computers from 20 years ago, while the language itself is more than half a century old. It has no way to make use of modern-day hardware features such as memory caching or concurrent programming. We did some direct comparisons (... a drag race) between the orbit propagator of GPSdancer (JAVA) against the equivalent propagator of ESOC's IGS analysis centre (FORTRAN 90), and then spent a considerable time poking fun at our colleagues.

## Category 1: orbit propagation models

### 1.1 Numerical integrator

Numerical integrators tend to look complicated, but they all boil down to the same principles of extrapolation of polynomials through $$n$$ earlier points. In GPS Dancer the satellite orbits are propagated by a single integrator for the velocity (Adams method) in combination with a double integrator for the position (Gauss-Jackson method, a.k.a. "second sum" method).

The Adams integrator uses a quadrature expression in terms of the accelerations and velocities from $$n$$ previous points to predict the velocity in the new step. The Gauss-Jackson integrator uses a quadrature expression in terms of the accelerations and positions from the n previous points to predict the position in the new point. This allows evaluation of the acceleratioin in the new point. It then applies a corrector quadrature that uses $$n$$ position and acceleration values up to and including the predicted point. This improves the position vector in the new point.

Both integrators are currently used to order 8 (meaning that $$n$$ = 9), but all required quadrature coefficients are in fact computed during program initialization from their known recurrent relations. This means that the integrators may one day be redefined to use order 10 or more, without any need to reprogram the coefficient tables. It also means that there cannot be any typos in hardcoded quadrature coefficients, which is a very common problem in complex numerical integrators (...try typing a few dozen numbers of 14 digits each without a single mistake).

Start-up of all quadratures is done by iterated quadratures that express any of the 9 previous points in terms of the other 8 points. Coarse initialization (as required to start-up the start-up quadratures...) is done by first order Runge-Kutta back-propagation from the known starting vector. In other words, the given starting vector of any orbit is the first proper state vector of the orbit arc. Virtual poiints are created for the points that preceed the known vector, after which the nominal integration process can propagate forward. This approach allows easy restarting of the integrator for e.g. a smaller step size (eclipses, manoeuvres) because n-1 new quadrature points (at smaller step size) can be created from any known point, in such a way that the last known point before the change in step size now becomes the n-th point of the new quadrature buffer.

For more details on both the Adams and Gauss-Jackson integrator, see a full description in Journal of Astronautical Sciences, Vol. 52, No. 3, July-Sep 2004 pp 331-357.

### 1.2 Earth gravity

GPS Dancer uses the EGM2008 model up to degree and order 12. For ease of portability and interfacing, the coefficients up to 12x12 are embedded (hardcoded) in the class Gravity.java. Note however that a JAVA run-time environment imposes a size limit of 64 kB on class sizes, so that expansion to much higher degree (e.g. 180x180 for LEO support) will require an interface to external coefficients from an input file.

There is no need to worry about the difference between the tide-free and zero-tide version. They are only different in the C20 coefficient, and the value from the coefficients file is overruled by the time-variable expression of table 6.2 in IERS Technical Note 36.

### 1.3 Solid Earth tides

Dancer uses the precise model from the IERS2010 standards, chapter 7 of TN 36. The model has been verified with a set of test values generated by the SOFA implementation, a practice that is used for almost all other models as well. The effects up to degree 4 are included by Dancer.

### 1.4 Permanent tide

Exactly the same comment applies as discussed above for the solid earth tides. The permanent tide and its relation to C20 and geocentre are explained much better in the 2010 conventions (TN 36) than in the 2003 standards - thank you to whoever sorted out that mess.

### 1.5 Ocean tides

Dancer uses the FES2004 ocean tide model. Effects up to degree 12 are included (because that is the maximum geopotential degree used for GNSS satellites).

Over short orbit propagation periods (e.g. 24 hours), switching the ocean tide model off hardly has any effect on the orbit. Presumably, if we do not use the model there will be some long-term effects showing up in the orbits, so we keep using it. However, there should hardly be any impact on the station polyhedrons produced by routine GPSdancer solutions if ocean tides were ignored completely.

### 1.6 Third body gravity

The effects of direct gravity from Sun, Moon and four relevant planets (Mars, Venus, Jupiter, Saturn) is modelled in exact agreement with the IERS 2010 standards. The other planets are either too small (Mercury) or too far away (Uranus, Neptune) to have any relevant effect on GPS orbits. Poor Pluto is too small and too far away, and is not even considered a planet anymore. For the planetary ephemerides, the DE421 polynomials are used, and tested against the testdata provided by JPL. Effects of barycentre motion are of course properly included in all computations of gravity and the celestial reference frame. Mass ratios are taken from the IERS 2010 standards (that is, compatible with the DE421 ratios).

Curiously, the IGS REPRO standards discuss albedo and Earth infra-red radiation pressure but not the much larger effects of solar radiation pressure on the GPS satellites. It is of course known which models are used by the various IGS analysis centres, and there are in fact remarkable differences. Dancer contains both the 6-parameter and 9-parameter variants of the empirical "Bernese" model. It is intended to experiment with both to see which provides the best results, knowing that the 6-parameter model is to be preferred on theoretical grounds.

### 1.8 Earth albedo and infra-red radiation pressure

These effects are quite new to GPS orbit modelling, but are important for scale of the reference frame. At the moment, Dancer does not yet include these effects but an update to full compatibility with the IGS REPRO 2 standards will mend that at some point in the future. It is not expected that there will be any noticable difference in individual solution arcs, but over long periods there will be an improvement in seasonal and annual spectral defects in the GPS orbit solutions.

### 1.9 Effects of general relativity on orbital dynamics

The model of the IERS 2010 conventions is followed, but for the TT transformation (from numerical integer count of seconds) the IERS 2003 standards are still followed This is because the new description in the IERS TN 36 seems to indicate that nothing has changed since the 2003 standards, and then comes with a new equation.

Some clarification is sought of this issue, but presumably any future change in the definition will be virtually invisible in the results (we have a < 1cm orbit fit to the IGS orbits, so GPSdancer cannot be so wrong).

### 1.10 Eclipse modelling

The eclipse model in Dancer is a native implementation, and purely based on the projection of the disks of the Moon and Earth relative to the disk of the Sun (all seen from the satellite position). The computation is relatively simple as long as effects of Earth flattening are ignored (they are, but they shouldn't...). In that case, the edges of the (circular) disks of Sun, Moon and Earth each define two crossing points along the Sun's orbit, as seen from the satellite, so that only these six points need to be computed. Eclipses occur as soon as any of the four points defined for Earth and Moon fall between the two points of the sun. Then some additional computations are done for computing the exact umbra/penumbra situation.

The power transmitted by the GPS payload is always aimed at the Earth, i.e. in the nadir direction of the satellite. This is a very systematic effect, and leads to minor error in the scale if not included in the dynamic models. GPSdancer does not yet model this, but an update is foreseen (the model is closely related to the albedo / IR radiation pressure discussed above).

### 1.12 Effects of thrusting (orbit manoeuvres)

Manoeuvres can be cleanly included in the GPSdancer orbit propagation without secondary numerical defects induced by discontinuities in higher derivatives of the acceleration function. This is accomplished by modelling a manoeuvre as a cosine-shaped acceleration, starting and ending at zero and with derivatives of zero. The integration step size is decreased substantially during pulse integration, and returns to nominal (in phase with the earlier steps) afterwards. Two empirical scale factors are estimated per manoeuvre, representing directly the delta-V in in-plane and out-of-plane directions.

Manoeuvre detection is not yet implemented in GPSdancer, but it will be. This is extremely important for near real-time processing, because failure to detect a satellite manoeuvre will easily lead to poor quality analysis results. Manoeuvre detection can be done effectively on the basis of observing a ramp function in the range residuals, i.e. fitting a straight line through the most recent 20 minutes of data and testing if the first derivative is above a certain threshold, while the overall noise (after removal of the ramp) is nominal.

## Category 2: GPS observation models

### 2.1 Station coordinate definition

As long as the receiver positions are estimated in the form of a displacement relative to the a priori coordinates, there is some freedom in defining the exact station marker position, although the obvious goal is to be strictly compliant with the common definitions as followed by IGS and IAU. However, all GPS Dancer processing is currently based on RINEX input data. For non-ITRF stations, the a priori coordinates are obtained from this source, even though the header field is optional and some RINEX files do not use it. This implies that GPSdancer will not work for those stations.

It is also possible to provide starting values in user configuration, but for automated processing this is not a good option (velocities tend to be inaccurate so that the values will gradually drift). In RINEX the definition of the a priori marker coordinates is in fact rather convoluted, because alternative approaches are still allowed in the header records, using offsets (delta) in either XYZ or East/North/Up directions, and at the same time binding the antenna phase centre to the clock via the absolute ANTEX antenna maps.

To clarify various misconceptions that exist regarding RINEX station coordinate definition, Dancer uses the a priori values of the APPROX POSITION XYZ header record, and applies DELTA  H/E/N to this vector if there is any such offset specified. The resulting point is assumed to be the antenna reference point, which is the origin of the antenna phase centre variations defined in ANTEX files. In other words, if your marker is not the antenna bore point (and usually it is not), you are encouraged to provide the offset in the DELTA H/E/N field, obviously as accurate as feasible. The indication "APPROX" is confusing: if the value is not accurate enough, Dancer may reject most of your data during initial screening.

### 2.2 Station displacements due to body tides

Dancer uses the exact body tide uplift model as specified in the IERS 2010 standards, chapter 7.1.

Dancer follows the FES2004 definition for ocean loading, but implements the version that requires the loading coefficients as external input. The available FORTRAN software for the complete model - allowing evaluation at arbitrary points on the planet - will at some point be ported to JAVA, but the hope exists that the people who have produced this terrible piece of code will at some point be so ashamed of their work that they produce a cleaned up version. We will then port the cleaned up version to JAVA. Until that time, you need to produce input coefficients for your receiver in the "BLQ" format, for which a website exists here. Please tick "YES" for the correction of the Earth's centre-of-mass.

### 2.4 Permanent tide

GPSdancer uses the IERS 2010 convention on the permanent tidal uplift. The entire issue of permanent tide, zero-tide and tide-free results is explained much better in the TN 36 than in previous standards. The modelled station coordinates are the zero-tide versions, because the tide-free coordinates are not observable. Therefore, the permanent tide is implicitly part of the coordinate definition. If we remove it, we would get tide-free coordinates, which have no practical meaning.

No effects of atmospheric loading are implemented in GPS Dancer (yet). Over the short 48-hour periods of Dancer analysis arcs this has no effect, but long time series of dancer station coordinates may reflect the absence of an atmospheric loading model. As with other models this will be fixed at some unknown point in the future. However, the IGS REPRO 2 standards also insist that no atmospheric loading models are included.

### 2.6 Absolute antenna phase centre

Dancer uses ANTEX antenna phase centre maps as inputs, and intends to always have the most recent version of the file. This file will be distributed by the Dancer peer-to-peer network itself (after being injected by a system administrator in at least one on-line peer) via automatic update facilities of the GPS Dancer software. As a user, you only need to specify the directory in which the ANTEX file resides, and of course which combination of receiver, antenna and dome you are using. If you use the Graphical User Interface, you will find a pull-down menu with all available options, based on the lists of formal receiver indications that is being maintained by the IGS.

For the satellites, the ANTEX maps are also used but no azimuth-dependence is implemented. If future ANTEX data defines azimuth dependence in addition to elevation dependence, the software will be modified accordingly.

### 2.7 Troposphere delay

The troposphere model in GPSdancer is now in complete agreement with the IERS 2010 conventions, using the vienese model for a priori dry and wet terms, and the associated mapping function in combination with estimated zenith delays and horizontal gradient parameters. Zenith delays are estimated once per hour, but the gradient parameters are estimated as a linear drift between a pair of parameter values (Teast, Tnorth), at the start and the end of the arc (i.e. four gradient parameters per arc).

The self-test option of GPSdancer runs a unit test on the test-values that are provided in the vienese (...FORTRAN) source code, and GPSdancer matches these results within 0.1% or so. Close enough.

### 2.8 Travel time effects

The 3D displacement of the satellite during the signal travel time is taken into account linearly, i.e. the spacecraft velocity vector is used to iteratively locate the satellite position at time of transmission. This is done after applying all other geometric corrections, including troposphere delay and statioin clock delay, but not the satellite clock delay. The station clock delay causes a small shift of the station in longitude direction, but also causes a minor displacement of the satellite.

### 2.9 Phase wind-up effects

For the effect of phase wind-up during a satellite's pass over a station, the precise model is used as described by Kouba in the IGS processing standards, section 5.1. For the satellite azimuth reference, the local North direction is used. The phase wind-up model is a notorious source of bugs in GPS analysis software, notably because its offset typically needs to be reset before reprocessing a given GPS pass. It also imposes the condition that all observations of a single station-satellite pass are processed in chronological order, which is not necessary for any of the other models.

### 2.10 Relativistic effects on signal propagation

Dancer follows the IERS 2010 standards for effects of relativity on signal propagation (chapter 11), but without worrying about the differences between barycentric and geocentric time, and without the gravitational term. The effects of this simplification were found to be smaller than 0.1 mm on GPS range, which is considered acceptable for the foreseeable future of GPS analysis.

### 2.11 Biases, biases, biases...

The plague of recent years in precise GNSS analysis is the arrival of a multitude of biases, that we could safely ignore ten or twenty years ago because all of our analysis was inaccurate anyway. Basically, there is a non-zero instrument bias between any two observation types, or between the same observation types at different carrier frequencies. These biases are not constant, but vary with time, receiver, satellite and antenna. If we are honest, they also vary from one analysis centre to another.

This is bad enough if you only process ionosphere-free combinations of code and phase (like GPS Dancer), which involves three biases, but two of them conveniently disappear (phase relative to reference code signal, on two frequencies) by estimating the phase ambiguity parameters. The remaining bias (difference between code signals on two carrier wave) has a receiver-dependent part and a transmitter-dependent part. The situation gets significantly worse with the arrival of a third frequency. Receiver biases are dependent on the brand and manufacture of your receiver, and satellite biases depend on particular block versions of the GPS satellite.

There is no absolute truth left in GPS analysis: a range measurement can be anything you like it to be, as long as you define an appropriate bias for it to fit observed reality. The current situation regarding bias definitions and computation in IGS is far from satisfactory.  An IGS Working Group has been created to do something about the situation, but this will probably take a considerable amount of time.

For the time being GPS Dancer does what everybody else does. We download a monthly file with differential code biases P1-C1 and C2-C1 from an obscure ftp server that is allegedly located in Berne, Switzerland. The file was put their after some non-documented magic mystery analysis by the CODE analysis centre, and we all trust that they follow some sort of intelligent process to arrive at the published numbers.

## Category 3: models used both in orbit propagation and observation modelling

### 3.1 Earth rotation model

The terrestrial reference frame realization in Dancer follows the IERS 2010 standards in every detail, i.e. also including all sub-daily tidal effects on polar motion, and what not. A set of test matrices was generated with trusted (IGS) software, and the GPSdancer Earth rotation model matches the external reference within 0.1 mm at the Earth's surface. This level of accuracy is actually necessary, if we aim for computing ITRF station coordinate time series with a noise level below 3 mm in 3D, and better than 1 mm horizontal.

Some comments need to be made regarding the definition of the Earth Rotation parameters in GPSdancer, given that nominal solution arcs span 48 hours and shift by e.g. 30 minutes relative to the previous solution. Dancer will always define a set of three Xp, Yp and dUT1 parameters at the start, middle and end of the 48-hour period. These estimated parameters are therefore not strictly defined at midnight UTC, as IGS and IERS want. The dUT1 offset at the start of the arc is of course fixed to its a priori value - absolute UT1 is unobservable from GPS, and Dancer can not yet handle VLBI data - so that there are always 8 estimated Earth rotation parameters for each 48-hour arc. Their values are initialized from the previous solution arc, accounting for linear drifts over the shift in arc epoch.

This approach is perfectly equivalent to IGS practice, as long as the Dancer outputs are not used in the IGS combination. If one day we want to do that, new parameters can be derived from the Dancer time series directly - i.e., we determine values at midnight crossings by some polynomial fit through the GPSdancer values.

### 3.2 GPS Satellite attitude model

The current implementation in Dancer matches more or less the standards from 10 years ago, but a lot has happened in the area of attitude modelling of GNSS satellites. In particular, behaviour of various individual satellites needs to be modelled (i.e. automatic attitude adjustents in eclipses, or maintaining sun-pointing mode for others) and the consistency of attitude modelling with e.g. antenna phase centre computations and phase wind-up can be improved. Several of these effects reache amplitudes of a few cm, so that they are not negligible.

As in all models, the aim is to comply with the REPRO 2 standards of IGS. In fact, we would encourage the IGS to consider the production of attitude quaternion files as a new combination product, which could really help many precise users to improve their results. Attitude quaternions are already very common in LEO missions such as GRACE, Jason, Sentinel, and make it much easier to ensure correct modelling of the antenna offsets at the spacecraft side.

### Project details

The GPS Dancer project started in 2007 as a voluntary project of a working group of the International Association of Geodesy.

### Square dance algorithm

The GPS Dancer system was named after its "square dance" exchange algorithm. Of course, it also wants to to make the GPS reference frame denser.