Numerical Error in Interplanetary Orbit Determination Software

The most used method for representing real numbers in modern computers is the binary floating point representation, based upon the normalized scientific notation in base $2$. Using this representation an exact value $x_{0}$ is encoded using a finite number of bits, as:

$$x=(-1)^{s}(1.b_{1}b_{2}\ldots b_{t-1}b_{t})\times 2^{p}$$

(1)

where $s=0$ for positive values, $s=1$ for negative ones. $p$ can be computed from $x_{0}$ using the following relation:

$$p(x_{0})=\begin{cases}\lfloor\log_{2}{|x_{0}|}\rfloor,&\text{if $\quad x_{0}\neq 0$}\\ 0,&\text{if $\quad x_{0}=0$}\end{cases}$$

(2)

where $\lfloor a\rfloor$ is the floor function of $a$.

$b_{t}$ is the last represented digit of $x$ and is called Least-Significant Bit ($LSB$). In the representation of $x_{0}$ the maximum value represented by the $LSB$ is:

$$q(x_{0})=2^{p(x_{0})-t}$$

(3)

There are several rounding algorithms that define how to choose $b_{t}$. The mostly used algorithm is “round toward nearest, ties to even”, because it minimizes the rounding error for a given $t$ and does not introduce a systematic drift Goldberg:1991 . In what follows, only this algorithm will be considered.

The use of a finite number of digits to represent the mantissa introduces an error. The absolute and relative rounding errors are defined as:

	$$\displaystyle\Delta\overline{x}(x_{0})\triangleq x-x_{0}$$		(4)
	$$\displaystyle\Delta\tilde{x}(x_{0})\triangleq(x-x_{0})/x_{0}$$		(5)

The maximum absolute value of the absolute error in the representation of a variable $x_{0}$ is:

$$\Delta\overline{x}_{\text{max}}(x_{0})\triangleq\max{|\Delta\overline{x}(x_{0})|}=q(x_{0})/2=2^{p(x_{0})-t-1}$$

(6)

Hence, the maximum absolute error depends upon the number of digits available for the fractional part of the mantissa $t$ and the order of magnitude $p$ of the binary number. Given $t$, the maximum error is the same for all values in the range $[2^{p},2^{p+1})$.

Given two proportional values $x_{1}$ and $x_{2}=k\,x_{1}$, in general $\Delta\overline{x}_{\text{max}}(x_{2})$ and $\Delta\overline{x}_{\text{max}}(x_{1})$ do not satisfy the same relation of proportionality, because of the presence of the floor function in the definition of $p$ (Eq. 2). If and only if the absolute value of the scale factor $k$ is an integer power of the basis of representation ($|k|=2^{n_{k}}$, $n_{k}\in\mathbb{Z}$) the maximum absolute errors are proportional, with scale factor $|k|$:

$$p(x_{2})=\lfloor\log_{2}{|kx_{1}|}\rfloor=\lfloor\log_{2}{|k|}+\log_{2}|x_{1}|\rfloor=\lfloor n_{k}+\log_{2}|x_{1}|\rfloor=n_{k}+\lfloor\log_{2}|x_{1}|\rfloor=n_{k}+p(x_{1})$$

(7)

$$\Delta\overline{x}_{\text{max}}(x_{2})=2^{p(x_{2})-t-1}=2^{n_{k}+p(x_{1})-t-1}=2^{n_{k}}\,2^{p(x_{1})-t-1}=|k|\,\Delta\overline{x}_{\text{max}}(x_{1})$$

(8)

However, it can be shown that for a generic scale factor $k$ the ratio $\Delta\overline{x}_{\text{max}}(x_{2})/\Delta\overline{x}_{\text{max}}(x_{1})$ is still $|k|$, when averaged on a large enough interval.

As an immediate consequence, the maximum absolute error in the floating point representation of a specific physical quantity depends upon the adopted unit. For example:

	$$\displaystyle\text{LSB}(1\,\text{km})\simeq 2.22\times 10^{-16}\,\text{km}=2.22\times 10^{-13}\,\text{m}$$		(9)
	$$\displaystyle\text{LSB}(1000\,\text{m})\simeq 1.14\times 10^{-13}\,\text{m}$$		(10)

However, averaging on a large enough interval of values, the two maximum absolute errors are equal, so they can be considered equivalent.

The maximum absolute value of the relative error, also called machine epsilon or $\varepsilon$, is:

$$\varepsilon\triangleq\max{\left|\Delta\tilde{x}(x_{0})\right|}=\Delta\overline{x}_{\text{max}}/|x_{\text{min}}|=2^{-t-1}$$

(11)

$\varepsilon$ does not depend upon the order of magnitude of the value $x_{0}$, but only upon the number of digits of the fractional part of the mantissa $t$.

The IEEE-754 standard defines several basic floating-point binary formats that differ by the number of bits used to encode the sign, the exponent and the fractional part of the mantissa. The most important formats are Binary32 (also called single precision), Binary64 (also called double-precision), and Binary128 (also called quadruple-precision). Their characteristics are described in Table 1 IEEE-754:2008 .

The rounding error is a form of quantization error, with the LSB as the quantization step. Hence, it is possible to adopt a statistical model, usually applied to quantization errors, known as uniform quantization error model, or Widrow’s model Widrow:1961 . Using this description, the round-off error $\Delta\overline{x}(t)$ in the numerical representation of an exact time function $x_{0}(t)$ is modeled as a white stochastic process, wide-sense stationary, uncorrelated with the represented function $x_{0}(t)$, and with a Probability Density Function (PDF) uniform between the extreme values $-q(x_{0})/2$ and $q(x_{0})/2$.

With this model, it is possible to compute the statistical characteristics of the numerical error:

1. 1)

   Mean value: $\operatorname{E}[\Delta\overline{x}(t)]=0$

2. 2)

   Maximum absolute value: $\Delta\overline{x}_{\text{max}}=q/2$

3. 3)

   Variance: $\sigma_{\Delta\overline{x}}^{2}=q^{2}/12$

The model is valid if the PDF and the Characteristic Function (CF) (the Fourier transform of the PDF) of the input signal $x_{0}(t)$ satisfy the conditions of the Quantization Theorem (QT), stated in Widrow:1961 . In Sripad:1977 weaker necessary and sufficient conditions for the QT are provided but, in practice, for real-world input signals, the conditions are very restrictive. For example, a sinusoidal signal does not satisfy the QT, but Widrow’s model can still be applied with sufficient accuracy if the amplitude of the signal is several times larger than the quantization step Mariano:2006 .

In this work the statistical model was assumed to be valid and this assumption was verified through a detailed analysis of the rounding errors, described in Sec. III.3.

When performing an exact operation $f(x_{0},y_{0})$ using floating point arithmetic, the result is affected by two kind of errors with respect to the theoretical value $z_{0}$:

1. 1)

   Propagation error of the rounding errors in the inputs:

   $$x=x_{0}+\Delta\overline{x},\qquad y=y_{0}+\Delta\overline{y}$$

   (12)

   If the errors are small enough:

   $$\tilde{z}_{0}=f(x,y)\simeq f(x_{0},y_{0})+\frac{\partial f}{\partial x}(x_{0},y_{0})\Delta\overline{x}+\frac{\partial f}{\partial y}(x_{0},y_{0})\Delta\overline{y}$$

   (13)

2. 2)

   Additional rounding error: even if the errors in the inputs are zero, the output could be affected by a rounding error. IEEE Standard 754 requires that the result of the basic operations (addition, subtraction, multiplication, division and square root) must be exactly rounded IEEE-754:2008 , i.e. the operation shall be performed as if it first produced an intermediate result correct to infinite precision, and then that result was rounded to the finite number of digits in use. The rounding on the result $\tilde{z}_{0}$ is represented by the addition of the random variable $\Delta\overline{z}$:

   $$z=\tilde{z}_{0}+\Delta\overline{z}=f(x,y)+\Delta\overline{z}$$

   (14)

   In practice, to obtain an exactly rounded result, it is necessary to perform the basic operations using at least 3 additional binary digits Goldberg:1991 . For non-basic operations, the error on the results could be larger and depends on how the operations are implemented. The additional rounding error may be zero if the intermediate result of the operation does not need rounding.

The total error in the result is the sum of the propagation error and the additional rounding error:

$$\Delta z=\frac{\partial f}{\partial x}(x_{0},y_{0})\Delta\overline{x}+\frac{\partial f}{\partial y}(x_{0},y_{0})\Delta\overline{y}+\Delta\overline{z}$$

(15)

Using Widrow’s model for numerical errors $\Delta\overline{x}$, $\Delta\overline{y}$, and $\Delta\overline{z}$, the total error $\Delta z$ is the sum of three random variables and its statistical characteristics derive from the characteristics of each round-off error.

According to Moyer’s DRD formulation, the unramped two-way or three-way Doppler observable is proportional to the mean time variation of the round-trip light-time (i.e. the mean range-rate divided by $c$), during the count time interval centered in the observable’s time-tag Moyer:2000 :

$$F_{2,3}(TT)=M_{2}f_{T}\left[\rho(TT+T_{c}/2)-\rho(TT-T_{c}/2)\right]/T_{c}$$

(16)

where $M_{2}$ is the spacecraft transponder turnaround ratio, which is the ratio of the transmitted frequency to the received frequency at the spacecraft. $\rho$ is the time interval between the reception time, referred to the electronics of the receiving ground station and expressed in Station Time³³3Station Time is a realization of the proper time at the Earth ground station. (ST), and the corresponding transmission time, at the electronics of the transmitting ground station, in ST:

$$\rho(t_{3}(\text{ST})_{R})=t_{3}(\text{ST})_{R}-t_{1}(\text{ST})_{T}$$

(17)

Neglecting second order terms, such as relativistic light-time delays, time scale transformations, electronic delays, and transmission media delays, the round-trip light-time can be written as:

$$\rho(t_{3}(\text{ST})_{R})\simeq r_{23}/c+r_{12}/c=1/c\cdot\left[\left(\mathbf{r_{3}}-\mathbf{r_{2}}\right)\cdot\left(\mathbf{r_{3}}-\mathbf{r_{2}}\right)\right]^{1/2}+1/c\cdot\left[\left(\mathbf{r_{2}}-\mathbf{r_{1}}\right)\cdot\left(\mathbf{r_{2}}-\mathbf{r_{1}}\right)\right]^{1/2}$$

(18)

where all state vectors are expressed in the Solar System barycentric space-time Frame Of Reference (FOR).

In Eq. 18 the position vectors $\mathbf{r_{3}}$ and $\mathbf{r_{1}}$ are computed from the planetary and lunar ephemeris and from the Earth-centered position of the ground station antenna:

	$$\displaystyle\mathbf{r_{3}}=\mathbf{r^{C}_{A3}}(t_{3})=\mathbf{r^{C}_{B}}(t_{3})-\mathbf{r^{E}_{M}}(t_{3})/(1+\mu_{E}/\mu_{M})+\mathbf{r^{E}_{A3}}(t_{3})$$		(19)
	$$\displaystyle\mathbf{r_{1}}=\mathbf{r^{C}_{A1}}(t_{1})=\mathbf{r^{C}_{B}}(t_{1})-\mathbf{r^{E}_{M}}(t_{1})/(1+\mu_{E}/\mu_{M})+\mathbf{r^{E}_{A1}}(t_{1})$$		(20)

Superscripts and subscripts designating location are explained in Table 2.

$\mathbf{r_{2}}$ is computed from the spacecraft ephemeris, obtained by integrating the equations of motion with respect to the Center Of Integration (COI), and from the planetary ephemeris:

$$\mathbf{r_{2}}=\mathbf{r^{C}_{P}}(t_{2})=\mathbf{r^{C}_{O}}(t_{2})+\mathbf{r^{O}_{P}}(t_{2})$$

(21)

The reception time in Ephemeris Time⁴⁴4Ephemeris Time is the coordinate time of the Solar System barycentric space-time frame of reference used in the adopted celestial ephemeris. ($T_{eph}$) Standish:1998 , $t_{3}$, is computed from the reception time in ST, $t_{3}(\text{ST})$, through different time transformations while $t_{2}$ and $t_{1}$ are computed from $t_{3}$ solving for the down-link and the up-link light-time problems, respectively.

Using the DRD formulation, the Doppler observable is computed as the difference of two round-trip light-times (Eq. 16). Qualitatively, this formulation has a high sensitivity to round-off errors because, as it is well known, the difference between two large and nearly equal values substantially increases the relative error, causing a loss of significance. Hence, errors in the round-trip light-time have a large influence on the Doppler observables.

At first approximation, the errors $\Delta\rho$ in the round-trip light-time at the start and at the end of the count time can be considered independent. From Eq. 16, an order of magnitude estimation of the corresponding effect in the Doppler observable is:

$$\Delta F_{2,3}\simeq M_{2}f_{T}\sqrt{2}\Delta\rho/T_{c}$$

(22)

The main sources of numerical errors in the round-trip light-time are:

1. 1)

   Representation of times: the rounding errors in time epochs $t_{3}$, $t_{2}$, and $t_{1}$ propagate only indirectly into the round-trip light-time. An error $\Delta t_{k}$ in time epoch $t_{k}$ affects the computation of a position vector $\mathbf{r^{i}_{j}}(t_{k})$, through the corresponding velocity vector $\mathbf{\dot{r}^{i}_{j}}$:

   $$\Delta\mathbf{r^{i}_{j}}=\mathbf{\dot{r}^{i}_{j}}\,\Delta t_{k}$$

   (23)

   At first approximation, the total effect on $\rho$ is:

   $$\Delta\rho\simeq(\dot{r}_{12}+\dot{r}_{23})\Delta t_{k}/c$$

   (24)

   Where $\dot{r}_{ij}$ is the range-rate, i.e. the time variation of the range, between $i$ and $j$.

   Currently both the ODP and AMFIN express the time variable $t_{k}$ as a single double-precision value of the time elapsed since a reference epoch. In the ODP, time is measured in seconds past J2000 (January 1st 2000, 12:00) Moyer:2000 , while in AMFIN, time is measured in days past January 1st 2000, 00:00 AMFIN:ALG . Fig. 1 shows the maximum absolute round-off error in time variables, for both the double-precision time representations adopted in the ODP and AMFIN. When the time variable increases relative to the reference epoch, the maximum rounding error increases with a piecewise trend, because it doubles when the binary order of magnitude $p(t)$ increases by one.

   

   The difference of 12 hours in the reference epochs of the two programs is negligible. Hence, as discussed in Sec. II, because of the different measurement units, the time round-off errors in these two OD codes are not exactly the same, but they have the same average magnitude, on a large enough time interval. For the ODP-like time representation, between July 2008 and January 2017, the maximum round-off error is about $30\,\text{ns}$. From Eq. 24, considering a range-rate of $30\,\text{km/s}$, the effect of the time error on the two-way range has an order of magnitude of $2\,\text{mm}$. From Eq. 22, the approximate effect on the two-way range-rate is about $4\times 10^{-2}\,\text{mm/s}$, for a count time of $60\,\text{s}$.

2. 2)

   Representation of distances: the rounding errors in each component of the position vectors propagate into the precision round-trip light-time directly, because the position vectors are used in the computation of the Newtonian distances $r_{12}$ and $r_{23}$, but also indirectly, because the position vectors are used also in the computation of $t_{2}$ and $t_{1}$. However, the indirect effects are much smaller than the direct ones and can be neglected. Both the ODP and AMFIN express distances in km, so at $1\,\text{AU}$ the maximum round-off error in the two-way range is about $0.03\,\text{mm}$ and at $10\,\text{AU}$ it is about $0.24\,\text{mm}$. From Eq. 22, the approximate effect on the two-way range-rate is about $7\times 10^{-4}\,\text{mm/s}$ at $1\,\text{AU}$ and $6\times 10^{-3}\,\text{mm/s}$ at $10\,\text{AU}$, for a count time of $60\,\text{s}$.

3. 3)

   Additional rounding errors: the rounding errors in the result of each operation propagate into the round-trip light-time. In the computation of the round-trip light-time $\rho$, the following computational steps were considered:

   	$$\displaystyle\mathbf{r^{E}_{B}}(t_{i})=\mathbf{r^{E}_{M}}(t_{i})/(1+\mu_{E}/\mu_{M}),\qquad i=1,3$$		(25)
   	$$\displaystyle\mathbf{r^{C}_{E}}(t_{i})=\mathbf{r^{C}_{B}}(t_{i})-\mathbf{r^{E}_{B}}(t_{i}),\qquad i=1,3$$		(26)
   	$$\displaystyle\mathbf{r_{i}}=\mathbf{r^{C}_{Ai}}(t_{i})=\mathbf{r^{C}_{E}}(t_{i})+\mathbf{r^{E}_{Ai}}(t_{i}),\qquad i=1,3$$		(27)
   	$$\displaystyle\mathbf{r_{2}}=\mathbf{r^{C}_{P}}(t_{2})=\mathbf{r^{C}_{O}}(t_{2})+\mathbf{r^{O}_{P}}(t_{2})$$		(28)
   	$$\displaystyle\mathbf{r_{ij}}=\mathbf{r_{j}}-\mathbf{r_{i}}=\left[x_{ij}\quad y_{ij}\quad z_{ij}\right]^{\text{T}},\qquad(i,j)=(1,2),(2,3)$$		(29)
   	$$\displaystyle xx_{ij}=x_{ij}^{2},\quad yy_{ij}=y_{ij}^{2},\quad zz_{ij}=z_{ij}^{2},\qquad(i,j)=(1,2),(2,3)$$		(30)
   	$$\displaystyle xy_{ij}=xx_{ij}+yy_{ij},\qquad(i,j)=(1,2),(2,3)$$		(31)
   	$$\displaystyle rr_{ij}=xy_{ij}+zz_{ij},\qquad(i,j)=(1,2),(2,3)$$		(32)
   	$$\displaystyle r_{ij}=\sqrt{rr_{ij}},\qquad(i,j)=(1,2),(2,3)$$		(33)
   	$$\displaystyle t_{ij}=r_{ij}/c,\qquad(i,j)=(1,2),(2,3)$$		(34)
   	$$\displaystyle\rho=t_{12}+t_{23}$$		(35)

   Moreover, the position vectors $\mathbf{r^{i}_{j}}$ used in Eqs. 25, 26, 27, and 28 are computed using the astronomical ephemeris (given as an input to the program), the integration of the equations of motion and the Earth-fixed location of the ground station antenna. For simplicity, the numerical errors introduced in these computational steps are neglected. This assumption was verified a posteriori through the validation of the numerical errors model.

As discussed in Sec. II.3, a numerical error $\Delta\overline{x}$ in a generic variable $x$ propagates into the output variable $z$ through the partial derivative $\partial z/\partial x$. Hence, the numerical noise in the round-trip light-time $\rho$ at reception time $t_{3}$ can be computed through the partial derivative of $\rho$ with respect to each considered error source. Neglecting second order terms, the resulting expression for the total numerical error on $\rho$ is:

$$\Delta\rho(t_{3})=(\dot{r}_{12}+\dot{r}_{23})\Delta\overline{t_{3}}/c+(\dot{r}_{12}-\dot{p}_{23})\Delta\overline{t_{2}}/c-\dot{p}_{12}\Delta\overline{t_{1}}/c\\ +(\mathbf{\hat{r}_{12}}-\mathbf{\hat{r}_{23}})/c\cdot\left[\Delta\overline{\mathbf{r^{O}_{P}}}(t_{2})+\Delta\overline{\mathbf{r^{C}_{O}}}(t_{2})+\Delta\overline{\mathbf{r^{C}_{P}}}(t_{2})\right]\\ +\frac{\mathbf{\hat{r}_{23}}}{c}\cdot\left[\Delta\overline{\mathbf{r^{E}_{A3}}}(t_{3})+\Delta\overline{\mathbf{r^{C}_{B}}}(t_{3})-\frac{\Delta\overline{\mathbf{r^{E}_{M}}}(t_{3})}{1+\mu_{E}/\mu_{M}}+\Delta\overline{\mathbf{r^{C}_{A3}}}(t_{3})+\Delta\overline{\mathbf{r^{C}_{E}}}(t_{3})-\Delta\overline{\mathbf{r^{E}_{B}}}(t_{3})+\Delta\overline{\mathbf{r_{23}}}\right]\\ -\frac{\mathbf{\hat{r}_{12}}}{c}\cdot\left[\Delta\overline{\mathbf{r^{E}_{A3}}}(t_{1})+\Delta\overline{\mathbf{r^{C}_{B}}}(t_{1})-\frac{\Delta\overline{\mathbf{r^{E}_{M}}}(t_{1})}{1+\mu_{E}/\mu_{M}}+\Delta\overline{\mathbf{r^{C}_{A1}}}(t_{1})+\Delta\overline{\mathbf{r^{C}_{E}}}(t_{1})-\Delta\overline{\mathbf{r^{E}_{B}}}(t_{1})-\Delta\overline{\mathbf{r_{12}}}\right]\\ +\frac{\Delta\overline{xx_{12}}+\Delta\overline{yy_{12}}+\Delta\overline{zz_{12}}+\Delta\overline{xy_{12}}+\Delta\overline{rr_{12}}}{2\,c\,r_{12}}+\frac{\Delta\overline{xx_{23}}+\Delta\overline{yy_{23}}+\Delta\overline{zz_{23}}+\Delta\overline{xy_{23}}+\Delta\overline{rr_{23}}}{2\,c\,r_{23}}\\ +\left(\Delta\overline{r_{12}}+\Delta\overline{r_{23}}\right)/c+\Delta\overline{t_{12}}+\Delta\overline{t_{23}}+\Delta\overline{\rho}$$

(36)

where $\dot{p}_{ij}$ is defined as:

$$\dot{p}_{ij}\triangleq\mathbf{\dot{r}_{i}}\cdot\mathbf{\hat{r}_{ij}}\qquad(i,j)=(1,2),(2,3)\\$$

(37)

The numerical error in two- or three-way Doppler observables with time-tag $TT$ can be expressed as:

$$\Delta F_{2,3}(TT)=M_{2}f_{T}\left[\Delta\rho(TT+T_{c}/2)-\Delta\rho(TT-T_{c}/2)\right]/T_{c}$$

(38)

Adopting the statistical model for each numerical error in Eq. 36, it is possible to evaluate the statistical properties of the total numerical error in the round-trip light-time $\rho$ and Doppler observable $F_{2,3}$. It follows that $\Delta\rho(t_{3})$ is a stochastic process with the following characteristics:

1. 1)

   White.

2. 2)

   Non stationary, because the coefficients multiplying every single numerical error are a function of time. However their time variation is typically very slow and, for a typical duration of a tracking pass, the process can be considered stationary.

3. 3)

   Gaussian-like distribution/PDF: from the central limit theorem the sum of $N$ independent and identically distributed random variables converges to a normally distributed random variable as $N$ approaches infinity. Because the different numerical errors are a finite number and not identically distributed, the expected PDF is only qualitatively bell-shaped.

4. 4)

   Zero mean.

5. 5)

   Variance: the total variance is the sum of each term’s variance:

   $$\sigma^{2}_{\Delta\rho}(t_{3})=\left(\frac{\dot{r}_{12}+\dot{r}_{23}}{c}\right)^{2}\sigma^{2}_{\Delta\overline{t_{3}}}+\left(\frac{\dot{r}_{12}-\dot{p}_{23}}{c}\right)^{2}\sigma^{2}_{\Delta\overline{t_{2}}}+\left(\frac{\dot{p}_{12}}{c}\right)^{2}\sigma^{2}_{\Delta\overline{t_{1}}}\\ +\left(\frac{x_{12}/r_{12}-x_{23}/r_{23}}{c}\right)^{2}\left[\sigma^{2}_{\Delta\overline{x^{O}_{P}}}(t_{2})+\sigma^{2}_{\Delta\overline{x^{C}_{O}}}(t_{2})+\sigma^{2}_{\Delta\overline{x^{C}_{P}}}(t_{2})\right]\\ +\left(\frac{y_{12}/r_{12}-y_{23}/r_{23}}{c}\right)^{2}\left[\sigma^{2}_{\Delta\overline{y^{O}_{P}}}(t_{2})+\sigma^{2}_{\Delta\overline{y^{C}_{O}}}(t_{2})+\sigma^{2}_{\Delta\overline{y^{C}_{P}}}(t_{2})\right]\\ +\left(\frac{z_{12}/r_{12}-z_{23}/r_{23}}{c}\right)^{2}\left[\sigma^{2}_{\Delta\overline{z^{O}_{P}}}(t_{2})+\sigma^{2}_{\Delta\overline{z^{C}_{O}}}(t_{2})+\sigma^{2}_{\Delta\overline{z^{C}_{P}}}(t_{2})\right]\\ +\left(\frac{x_{23}}{c\,r_{23}}\right)^{2}\left[\sigma^{2}_{\Delta\overline{x^{E}_{A3}}}(t_{3})+\sigma^{2}_{\Delta\overline{x^{C}_{B}}}(t_{3})+\frac{\sigma^{2}_{\Delta\overline{x^{E}_{M}}}(t_{3})}{(1+\mu_{E}/\mu_{M})^{2}}+\sigma^{2}_{\Delta\overline{x^{C}_{A3}}}(t_{3})+\sigma^{2}_{\Delta\overline{x^{C}_{E}}}(t_{3})+\sigma^{2}_{\Delta\overline{x^{E}_{B}}}(t_{3})+\sigma^{2}_{\Delta\overline{x_{23}}}\right]\\ +\left(\frac{y_{23}}{c\,r_{23}}\right)^{2}\left[\sigma^{2}_{\Delta\overline{y^{E}_{A3}}}(t_{3})+\sigma^{2}_{\Delta\overline{y^{C}_{B}}}(t_{3})+\frac{\sigma^{2}_{\Delta\overline{y^{E}_{M}}}(t_{3})}{(1+\mu_{E}/\mu_{M})^{2}}+\sigma^{2}_{\Delta\overline{y^{C}_{A3}}}(t_{3})+\sigma^{2}_{\Delta\overline{y^{C}_{E}}}(t_{3})+\sigma^{2}_{\Delta\overline{y^{E}_{B}}}(t_{3})+\sigma^{2}_{\Delta\overline{y_{23}}}\right]\\ +\left(\frac{z_{23}}{c\,r_{23}}\right)^{2}\left[\sigma^{2}_{\Delta\overline{z^{E}_{A3}}}(t_{3})+\sigma^{2}_{\Delta\overline{z^{C}_{B}}}(t_{3})+\frac{\sigma^{2}_{\Delta\overline{z^{E}_{M}}}(t_{3})}{(1+\mu_{E}/\mu_{M})^{2}}+\sigma^{2}_{\Delta\overline{z^{C}_{A3}}}(t_{3})+\sigma^{2}_{\Delta\overline{z^{C}_{E}}}(t_{3})+\sigma^{2}_{\Delta\overline{z^{E}_{B}}}(t_{3})+\sigma^{2}_{\Delta\overline{z_{23}}}\right]\\ +\left(\frac{x_{12}}{c\,r_{12}}\right)^{2}\left[\sigma^{2}_{\Delta\overline{x^{E}_{A3}}}(t_{1})+\sigma^{2}_{\Delta\overline{x^{C}_{B}}}(t_{1})+\frac{\sigma^{2}_{\Delta\overline{x^{E}_{M}}}(t_{1})}{(1+\mu_{E}/\mu_{M})^{2}}+\sigma^{2}_{\Delta\overline{x^{C}_{A1}}}(t_{1})+\sigma^{2}_{\Delta\overline{x^{C}_{E}}}(t_{1})+\sigma^{2}_{\Delta\overline{x^{E}_{B}}}(t_{1})+\sigma^{2}_{\Delta\overline{x_{12}}}\right]\\ +\left(\frac{y_{12}}{c\,r_{12}}\right)^{2}\left[\sigma^{2}_{\Delta\overline{y^{E}_{A3}}}(t_{1})+\sigma^{2}_{\Delta\overline{y^{C}_{B}}}(t_{1})+\frac{\sigma^{2}_{\Delta\overline{y^{E}_{M}}}(t_{1})}{(1+\mu_{E}/\mu_{M})^{2}}+\sigma^{2}_{\Delta\overline{y^{C}_{A1}}}(t_{1})+\sigma^{2}_{\Delta\overline{y^{C}_{E}}}(t_{1})+\sigma^{2}_{\Delta\overline{y^{E}_{B}}}(t_{1})+\sigma^{2}_{\Delta\overline{y_{12}}}\right]\\ +\left(\frac{z_{12}}{c\,r_{12}}\right)^{2}\left[\sigma^{2}_{\Delta\overline{z^{E}_{A3}}}(t_{1})+\sigma^{2}_{\Delta\overline{z^{C}_{B}}}(t_{1})+\frac{\sigma^{2}_{\Delta\overline{z^{E}_{M}}}(t_{1})}{(1+\mu_{E}/\mu_{M})^{2}}+\sigma^{2}_{\Delta\overline{z^{C}_{A1}}}(t_{1})+\sigma^{2}_{\Delta\overline{z^{C}_{E}}}(t_{1})+\sigma^{2}_{\Delta\overline{z^{E}_{B}}}(t_{1})+\sigma^{2}_{\Delta\overline{z_{12}}}\right]\\ +\frac{\sigma^{2}_{\Delta\overline{xx_{12}}}+\sigma^{2}_{\Delta\overline{yy_{12}}}+\sigma^{2}_{\Delta\overline{zz_{12}}}+\sigma^{2}_{\Delta\overline{xy_{12}}}+\sigma^{2}_{\Delta\overline{rr_{12}}}}{(2\,c\,r_{12})^{2}}+\frac{\sigma^{2}_{\Delta\overline{xx_{23}}}+\sigma^{2}_{\Delta\overline{yy_{23}}}+\sigma^{2}_{\Delta\overline{zz_{23}}}+\sigma^{2}_{\Delta\overline{xy_{23}}}+\sigma^{2}_{\Delta\overline{rr_{23}}}}{(2\,c\,r_{23})^{2}}\\ +\left(\sigma^{2}_{\Delta\overline{r_{12}}}+\sigma^{2}_{\Delta\overline{r_{23}}}\right)/c^{2}+\sigma^{2}_{\Delta\overline{t_{12}}}+\sigma^{2}_{\Delta\overline{t_{23}}}+\sigma^{2}_{\Delta\overline{\rho}}$$

   (39)

   Where the variance of each numerical error can be computed using the formulation presented in Sec. II.2.

The numerical error $\Delta F_{2,3}$ has similar characteristics to $\Delta\rho$ but it is not white, because two consecutive Doppler observables are a function of the round-trip light-time at their mid-time. In fact, the count times are consecutive, so the end of a count time is the start of the next count time:

$$TT_{k}+T_{c}/2=TT_{k+1}-T_{c}/2$$

(40)

Given Eq. 40, two consecutive Doppler observables can be written as:

	$\displaystyle F_{2,3}(TT_{k})$	$\displaystyle=M_{2}f_{T}\left[\rho(TT_{k}+T_{c}/2)-\rho(TT_{k}-T_{c}/2)\right]/T_{c}$
		$\displaystyle=M_{2}f_{T}\left[\rho(TT_{k+1}-T_{c}/2)-\rho(TT_{k}-T_{c}/2)\right]/T_{c}$		(41)
	$\displaystyle F_{2,3}(TT_{k+1})$	$\displaystyle=M_{2}f_{T}\left[\rho(TT_{k+1}+T_{c}/2)-\rho(TT_{k+1}-T_{c}/2)\right]/T_{c}$
		$\displaystyle=M_{2}f_{T}\left[\rho(TT_{k+2}-T_{c}/2)-\rho(TT_{k+1}-T_{c}/2)\right]/T_{c}$		(42)

Hence, the rounding error in $\rho(TT_{k+1}-T_{c}/2)$ affects both $F_{2,3}(TT_{k})$ and $F_{2,3}(TT_{k+1})$. The autocorrelation function of $\Delta F_{2,3}$ is:

$$R_{\Delta F}(\tau)\triangleq\frac{\operatorname{E}\left[\left(\Delta F(t)\right)\left(\Delta F(t+\tau)\right)\right]}{\sigma_{\Delta F(t)}\sigma_{\Delta F(t+\tau)}}\\ \simeq\frac{2R_{\Delta\rho}(\tau)-R_{\Delta\rho}(\tau-T_{c})-R_{\Delta\rho}(\tau+T_{c})}{2\left[1-R_{\Delta\rho}(T_{c})\right]}=\begin{cases}1,&\text{if $\tau=0$,}\\ -1/2,&\text{if $|\tau|=T_{c}$,}\\ 0,&\text{if $|\tau|>T_{c}$}.\end{cases}$$

(43)

The variance of $\Delta F_{2,3}$ is proportional to the variance of $\Delta\rho$:

$$\sigma^{2}_{\Delta F_{2,3}}(TT)=\left(M_{2}f_{T}/T_{c}\right)^{2}\left[\sigma^{2}_{\Delta\rho}(TT+T_{c}/2)+\sigma^{2}_{\Delta\rho}(TT-T_{c}/2)\right]\simeq 2\left(M_{2}f_{T}\sigma_{\Delta\rho}(TT)/T_{c}\right)^{2}$$

(44)

Thus, the standard deviation of $\Delta F_{2,3}$, representing the numerical noise level in Doppler observables, depends on:

1. 1)

   Floating point word’s length: the number of bits used for the fractional part of the mantissa defines the amplitude of all round-off errors. At present, almost all hardware and compilers implement the IEEE-754 standard.

2. 2)

   Time: the time at which the OD problem is settled defines the order of magnitude of all time epochs, thus defining the quantization step and the amplitude of the round-off errors in each time variable.

3. 3)

   Geometry: the state vectors of all participants to the OD problem define two critical factors: the order of magnitude of the three components of each position vector, from which the range and the additional round-off errors derive, and the velocities of the participants that define the influence of time round-off errors in the observables.

4. 4)

   Count time: $\Delta\rho$ does not depend upon $T_{c}$, hence $\Delta F_{2,3}$ (and its standard deviation) varies with $T_{c}^{-1}$. As a reference a typical white-phase reference noise varies with $T_{c}^{-1/2}$.