Relativistic positioning: including the influence of the gravitational action of the Sun and the Moon and the Earth's oblateness on Galileo satellites

Uncertainties in the satellite world lines lead to dominant positioning errors. In the present work, using the approach presented in \cite{neu14}, a new analysis of these errors is developed inside a great region surrounding Earth. This analysis is performed in the framework of the so-called Relativistic Positioning Systems (RPS). Schwarzschild metric is used to describe the satellite orbits corresponding to the Galileo Satellites Constellation. Those orbits are circular with the Earth as their centre. They are defined as the nominal orbits. The satellite orbits are not circular due to the perturbations they have and to achieve a more realistic description such perturbations need to be taken into account. In \cite{neu14} perturbations of the nominal orbits were statistically simulated. Using the formula from \cite{col10a} a user location is determined with the four satellites proper times that the user receives and with the satellite world lines. This formula can be used with any satellite description, although photons need to travel in a Minkowskian space-time. For our purposes, the computation of the photon geodesics in Minkowski space-time is sufficient as demonstrated in \cite{neu16}. The difference of the user position determined with the nominal and the perturbed satellite orbits is computed. This difference is defined as the U-error. Now we compute the perturbed orbits of the satellites considering a metric that takes into account the gravitational effects of the Earth, the Moon and the Sun and also the Earth oblateness...


Introduction
In this work, the numerical codes developed in Puchades and Sáez (2014) are used to compute the positioning errors due to uncertainties in the satellite world lines, which are defined as the U-errors. In Puchades and Sáez (2014) the U-errors were computed as the difference in a user positioning when describing the satellite world lines with Schwarzschild metric and a statistical perturbation of those world lines. Now the description of the perturbed satellite world lines is made using a metric that takes into account the gravitational effects of the Earth, the Moon, the Sun, as well as the Earth oblateness (Earth quadrupole effect). This analysis is performed in the framework of the so-called Relativistic Positioning Systems (RPS). Both the satellite world lines and the user location are numerically computed and then the corresponding U-errors are calculated. Part of the work shown here was presented in Fullana et al. (2019).
Firstly, a short introduction on the difference between traditional positioning with Global Navigation Satellite Systems (GNSS) and positioning with RPS is presented. Basically, for positioning, four satellites are needed to locate a user in spacetime; i.e. a satellite for each unknown. The system of algebraic equations that has to be solved to calculate the user location is given by Eqs. (7) (see section 4 for more detail). These equations are based on the null line element travelled by the photon from the emission to the reception event (see Juang and Tsai (2009)). The system given by Eqs. (7) has zero, one or two possible solutions. A study of such solutions in RPS can be seen in Puchades and Sáez (2012). In the present work we assume that only one solution is present.
The current GNSS are based on a Newtonian approach including relativistic corrections. The GNSS would work ideally if all the satellites and the receiver were at rest in an inertial reference frame. Depending on the desired accuracy for the positioning, the necessary corrective terms have to be included (see Gomboc et al. (2013)). However, RPS are based on a relativistic treatment from the beginning. For an account of this approach see Coll et al. (2006a) and Coll et al. (2006b). Important work has been done using different relativistic models to describe the motion of artificial satellites. Brumberg and Kopejkin (1989) constructed the harmonic, dynamically non-rotating reference system for any body of the solar system and derived the equations of motion of test particles in the vicinity of the given body using this reference system. This technique was applied to the Earth and its satellites. Meanwhile Damour et al. (1994) introduced a particular General Relativity (GR) celestial mechanics framework. They computed the equations needed for developing a complete relativistic theory of artificial Earth satellites. They claimed that their approach was more satisfactory than the previous ones especially with regard to its consistency, completeness and flexibility. In Kostić et al. (2015) a model of RPS is presented in a more realistic space-time near the Earth with all important gravitational effects: Earth multipoles up to 6th order, the Moon, the Sun, Jupiter, Venus, solid and ocean tides, and Kerr effect. A recent approach was proposed by Roh et al. (2016), Roh (2018). They implemented a full set of first-order Post-Newtonian (PN) corrections in the high-fidelity orbit propagation KASIOP (Korea Astronomy and Space Science Institute Orbit Propagator). And they numerically evaluated their effects on orbital elements for the laser geodynamics satellite (LAGEOS) and laser relativity satellite (LARES) orbits.
In the paper Philipp et al. (2018) the relativistic orbital effects were estimated, which need to be considered for GRACE satellite orbits. They compared the magnitude of relativistic corrections in a GR space-time with PN corrections with various non-gravitational perturbations of satellite orbits and their results proved that for a GRACE satellite orbit (low Earth orbit), the relativistic acceleration is of the same order of several environmental effects. Hence, relativistic effects need to be considered for high precision space missions. The altitude for GRACE is approx. 456 Km and for Galileo approx. 23222 Km, then the order of magnitude of all orbit effects for Galileo is slightly different from that of GRACE. a) Connection with previous works.
In this research, some results about RPS, by Coll et al. (2010), are applied in order to locate a user in the Earth vicinity up to 10 5 Km, using the numerical codes developed in Puchades and Sáez (2014).
In the Schwarzschild case, the time-like geodesic equations for a satellite were numerically solved in Delva and Olympio (2009). The authors presented a preliminary study for a satellite that moves in the equatorial plane. In this case, they calculated the emission coordinates, from the inertial ones, for a few users with the same spatial coordinates who receive the satellite signal at different times. The same problem was studied using two other methods to compare the accuracies and the computational time. Moreover, inCadez et al. (2010), methods were described, in the Schwarzschild space-time, to find the emission coordinates from the inertial ones and vice-versa. In Puchades and Sáez (2016), some aspects of the method described inCadez et al. (2010), to obtain the inertial coordinates (positioning) from the emission coordinates, were used; for example, a first-order approximation was used, in the parameter GM ⊕ /r, of the time transfer function (see Teyssandier and Le Poncin-Lafitte (2008)). Moreover, in Puchades and Sáez (2016), other aspects such as the use of the analytical solution of Coll et al. (2010) at zero-order calculations (which are very efficient) were considered. Also in Puchades and Sáez (2016), the positioning errors associated to the simplifying assumption that photons move in Minkowski space-time (S-errors) were estimated and it was shown that the approach based on the assumption that the Earth's gravitational field produces negligible effects on photons may be used in a large region surrounding Earth.
b) The present work. In this paper, the computation of a satellite world line is performed. This calculation is done in a GR space-time metric (see Resolutions of IAU (2000) for more detail). The metric taken into account considers the effect of the Earth, the Moon, the Sun and the Earth oblateness. Solar radiation pressure and other non-gravitational effects are not considered although some of them have important contributions (see for instance Roh et al. (2016), Roh (2018)). However, the order of magnitude of their contributions depends on the distance to the Earth. Notice that here linear perturbation theory is not considered, as Kostić et al. (2015) do. A metric is used from the first step to describe the space-time. Then satellite geodesic equations are computed. Afterwards a Runge-Kutta algorithm is implemented to solve the geodesics ordinary differential equations. Although we are aware that other effects should be considered we have decided to start with the three effects commented in this paper to implement our algorithm. Later on we will continue incorporating more PN contributions and other terms. Great accuracy in the digits is needed when considering all PN contributions. The terms considered in the present work are sufficient to test the algorithm and prepare it for our future purposes in RPS. Our results show that the greatest effect comes from the Earth quadrupole, the second one from Moon and the third one from Sun (at the height of Galileo and GPS satellites). Those satellite world lines are needed to position a user in RPS.
The perturbations computed here using metrics improve our previous works based on statistical methods as: 1) A better description of the real satellite world lines is achieved. 2)The effect of each perturbing contribution in the satellite world lines is studied. 3) Also the combination of two of the three terms in the metric is studied and the three of them together. So, the orbits of the satellites are described depending on the terms considered. 4) Therefore, the contribution of each effect on the user's positioning can also be studied. 5) The value of the U-errors is now smaller. 6) That means a more precise calculation of the user's positioning.
c) Schedule of this paper. Firstly, in section 2, the satellite world lines for a Galileo satellite (fixed orbital radius) are computed using the metrics of this paper. A study of the behaviour of such orbits is made. Afterwards, in section 3, the satellite world line (in a space-time given by a metric with the gravitational potential of the Sun or the Moon or with the Earth's oblateness or all effects together) is computed for a satellite at different orbital radius from the geocenter. In section 4, the satellite geodesics recently calculated are used to obtain the user location in space-time with RPS. The positioning errors are also calculated. Different users located on different spherical surfaces are considered. These spheres are centred in the geocenter and take different radii (from Earth radius to 5 × 10 4 Km). In Appendix A, a study of the order of magnitude of the orbit effects taken here is computed, that is to say, Moon and Sun gravitational potential and Earth oblateness. These are some of the greatest effects on the orbits. In Appendix B, a short description of circular orbits of Galileo satellites is written. In Appendix C, the time-like geodesic equations (calculating the Christoffel symbols) are written from the metric that describes a space-time, with Moon gravitational potential or Sun gravitational potential or Earth oblateness or taking into account all effects together. Finally, in section 5, the main conclusions and perspectives are commented. d) Notation and general considerations. Some aspects about the notation of this article are now commented. The Greek indices take values from 0 to 3 and the Latin indices go from 1 to 3 and x 0 ≡ ct, where t is the time coordinate and c the light velocity. Also, we use G for the Gravitational constant and M A for the mass of the body A. The symbols (subscripts) ⊕, ⊙ and are assigned to Earth, Sun and Moon, respectively. We use this notation throughout the article. The symbol AU is for one Astronomical Unit. The Minkowski tensor η µν is defined as: η µµ ≡ (−1, 1, 1, 1) for the diagonal elements and null for the rest of components. In our codes and in the formulae, we choose the units in such a way than the light velocity is c = 1.
Moreover, throughout this paper, the reference system considered is the Geocenter Reference System (GRS). An appropriate method is used to represent some quantities in 3D, t = constant (user time), space-time sections. Colour bars and an appropriate pixelization are necessary. In this paper, as in Puchades and Sáez (2012), Puchades and Sáez (2014) and Puchades and Sáez (2016), the HEALPIx (hierarchical equal area isolatitude pixelisation of the sphere) package (Górski et al. (1999)) is used, with the same parameters and projections than in the paper Puchades and Sáez (2014).
We simulate the world lines of the Galileo background configurations and we use these world lines as initial conditions for the ordinary differential equations (ODE). The Galileo constellation is composed by 27 satellites (n s = 27), located in three equally spaced orbital planes (9 uniformly distributed satellites in each plane). The inclination of these planes is α in = 56 deg and the altitude of the circular orbits is h = 23222 Km; thus, the orbital period is close to 14.2 h. The satellites are numerated in such a way that the satellites 1 to 9, 10 to 18, and 19 to 27 are placed in distinct consecutive orbital planes. Initially, the trajectories are assumed to be circumferences whose centres are located in the Earth centre, which is also the origin of the almost inertial reference system used for positioning. In sections 2 and 3, the Galileo satellite used is number 1.

Satellite world lines (Galileo Constellation)
This section shows the Galileo satellite world lines computation with our algorithm. The description of the satellite world lines with Schwarzschild metric gives circumferences centred in the centre of the Earth. The Earth is approximated as spherically symmetric and with no rotation. These trajectories are defined as nominal trajectories. A short description of the circular satellite world lines can be seen in Appendix B. In the present work, the Schwarzschild satellite world lines are perturbed to compute the effects of the Sun, the Moon and the Earth quadrupole. A detailed description of the solution of the Ordinary Differential Equations (ODE) to compute the satellite world lines in this metric can be seen in Appendix C. The numerical solution of such ODE is also explained. This section is divided in two subsections. The first one 2.1, presents computations taking into account only the Moon and Sun gravitational fields and their sum. The second one 2.2 presents the effects of the Earth quadrupole on the orbits too. This second subsection includes the computations of the effect on the satellite of the Earth oblateness alone and of this oblateness plus the sum of the effects of Sun and Earth. In all computations, the alignment of Moon and Sun respect to the satellite has been taken into account when computing the orbit of the satellite. Notice that the figures represent the deviation from the corresponding circular orbits, nominal orbits, that would follow the satellites if the Earth's field would describe them, Schwarzschild metric. In such figures, the circular orbits would be straight horizontal lines.
Let us remark that this section and the next one are introduced so as to clearly describe the perturbing effects that are considered in this research. Once those effects have been presented, we use them to define the U-errors in a similar way as it was performed in Puchades and Sáez (2014) but now with the new metrics considerations. That is to say, the perturbations are not computed with a statistical algorithm but with some physical effects that affect the satellites. This represents an advance from the previous work as it can describe more realistic cases.

Sun and Moon effects
A gravitational term (Sun or Moon) can be considered as a satellite orbit perturbation. The B gravitational term reads: where B is the celestial body (B = , ⊙), either Sun or Moon (in our case). The first spatial derivatives with respect to the additional B potential stand: The initial conditions for the Moon and Sun are given from the Ephemeris in the GRS (reference plane: equatorial and rectangular coordinates). The Miriade Ephemeris Generator is used. Therefore, the initial conditions (velocity and position vectors) for the Moon and the Sun are the corresponding Ephemeris for the time 2018 − 12 − 13T 17 : 00 : 00.00. See Resolutions of IAU (2000) for more details about initial conditions. Nevertheless, the initial conditions for the satellite (velocity and position vectors) considered here are a Galileo satellite moving in a circular orbit with angular, Ω, and linear, v, velocity from Ω = GM⊕ R 3 , v = GM⊕ R . The following factor is also considered: with: Notice that the relativistic factor Eq.(4), enhances the GR treatment with the PN approximation (see Eq. (A.36), p. 210 of Gourgoulhon (2013)): Fig. 1 plots the radial distance, R, of a Galileo satellite, for two orbital periods, from the geocenter. There are three different computations: a) Considering only the Moon effect (blue). b) Considering only the Sun effect (red). c) Considering the Moon plus the Sun effect together (magenta). Moon and Sun contributions are considered in the metric. Fig. 1 allows us to compute the variation of the distance of the satellite to the geocenter, this indicates the deviation from the circular orbit. The circular orbit is the one described with the metric considering only the Earth contribution: Schwarzschild metric. In the figure it would correspond to a straight horizontal line that would start from the point where the other lines start from left to right. Therefore the effect of each perturbation alone or combined on the satellite orbit can be measured. One way to compute this variation is to calculate the difference between two different positions of the satellite in its orbit. The maximum variation of the distance to the geocenter will be the difference between the farthest position and the nearest one. Notice that absolute values are computed. R max is this maximum distance (maximum in the plot) and R min is the minimum one (minimum in the plot) when two orbital periods are considered. So, the maximum variation of the radial distance of the satellite to the geocenter, when two orbital periods are considered is R max − R min .
Next the values of this variation, R max − R min , for the cases we have considered are here presented as it can be seen in Fig. 1: • Around 600 m when only the gravitational contribution from Moon is considered.
• Around 200 m when only the gravitational contribution from Sun is considered.
• Around 700 m when the gravitational contributions from Sun plus Moon are considered together.
The order of magnitude obtained is the one expected for every quantity computed in this paper as (see the book Teunissen and Montenbru (2015), Chapter 3, and Roh (2018)): • At several hundreds of meters the second-largest effect producing orbit changes is the acceleration caused by Sun plus Moon.
• In spite of the fact that Sun is much more massive, the tidal acceleration in the geocentric frame is bigger for the Moon. This fact is due to the much lower distance of Moon to the satellite (according to Table 1, column 4, Appendix A).

Adding the Earth oblateness
The gravitational potential for Earth oblateness is (see Montenbruck and Gill (2005) and Teunissen and Montenbruck (2015): where the Legendre polynomial of degree 2 is P 2 (cos θ) = (3 cos 2 θ − 1)/2. The value of the zonal gravitational coefficient of degree 2 of the Earth J 2 is J 2 = 1.08263 × 10 −3 and θ and R ⊕ are the co-latitude of the satellite and the equatorial radius of the Earth, respectively.
The first spatial derivatives due to Earth oblateness are: where A(i) is 1 for i = 1, 2 and is 3 for i = 3. Fig. 2 plots radial distance, as a function of proper time (for 2 orbital periods), of a Galileo satellite from the geocenter (see previous subsection for the definitions of the quantities used here). In this subsection, we only present the effects considered in the metric of i) the Earth oblateness and ii) Moon plus Sun plus Earth oblateness. The study of Fig. 2, gives the maximum variation of the radial distance of the satellite, R max − R min (as defined in the previous subsection), during two orbital periods: • About 2 Km for Earth oblateness.
• About 3 Km for Earth oblateness plus direct tides from Sun and Moon.
Again we obtain the expected order of magnitude for the effects considered in this paper: the greatest effect on GNSS satellites is mainly due to Earth's oblateness (see the book Teunissen and Montenbruck (2015) and Roh (2018)). That is, effects at kilometre level are obtained.
We now compare the results of Figs. 1 and 2. As it can be observed in Figs. 1 and 2, the recovering of the satellite radial distance in 1 or 2 orbital periods is achieved in the following cases: • Direct tides from Sun (see Fig. 1), • Earth oblateness (see Fig. 2), that is to say, when direct tides from Moon are not considered. Satellite position shifting is obtained when Moon is considered.
Notice that the order of magnitude of the R-axes is different in both figures. In Fig. 1 the difference between each division is of 100 meters while in Fig 2 each difference is of 500 meters. This is done so because the order of magnitudes of the effects are different, greater when the oblateness is introduced, as it can be seen in the plots. Also notice that the presence of the Moon causes the bump in the effects from Moon + Sun (magenta in Fig. 1). The reason for not observing this bump in Fig. 2 in the Moon + Sun + Earth oblateness (violet in Fig. 2, where Moon effect is also considered) is this difference in the order of magnitude of divisions of R-axes. We are more interested in emphasizing orders of magnitude and that's why this bump is smoothed.

Orbital perturbation effects for different orbital radius from the geocenter
In this section, the satellite world lines are calculated with our algorithm for different satellite orbital radii from the geocenter. So we extend the computation to other possible satellite orbits far from Galileo ones. The orbital perturbation effects are estimated for these satellite orbital radii. Again, the radial distance R versus proper time is evaluated for different satellite orbital perturbating effects: presence of Moon and Sun gravitational field and Earth quadrupole.

Considering different orbital perturbation effects for a given satellite altitude
In this subsection, all the orbital perturbation effects considered are evaluated for a given satellite altitude.
Figs. 3 and 4 show the radial distance of a satellite versus the proper time (for two orbital periods), from the geocenter, at satellite orbital radius 5 × 10 4 Km, 10 5 Km and 1.5 × 10 5 Km. The perturbing orbit effects that are considered in the metric are Earth oblateness, Moon gravitational potential, Sun gravitational potential, Moon plus Sun gravitational potential and, Earth oblateness plus Moon plus Sun gravitational potential.
From Fig. 3, we can conclude that when the satellite is orbiting at 5 × 10 4 Km the perturbing effects considered such as Earth oblateness, Moon gravitational potential or Sun gravitational potential are of the same order of magnitude, but the Moon's effect is the strongest one (see the amplitude between maximum and minimum distance to the geocenter for the cases when the Moon gravitational potential is included). Note that the Earth oblateness effect is bigger than Sun and Moon effect at the height of Galileo satellites, but this is not the case at 5 × 10 4 Km of orbital radius. All these results are compatible with the figure which gives the order of magnitude of various perturbations of a satellite orbit from the book Montenbruck and Gill (2005).
From Fig. 4, we can conclude that as the orbital radius of the satellite increases, the Earth oblateness effect decreases and the Moon and Sun effects increase. When the satellite is approaching the Moon, the Moon effect is bigger than the Sun one.

Considering one orbital perturbation effect varying the satellite altitude
In this subsection, different orbital perturbation effects (Moon plus Sun gravitational field and Earth quadrupole) are evaluated for different satellite altitudes. Fig. 5 shows the radial distance, divided by its given orbital radius, for a satellite, versus the proper time (for two orbital periods), at different orbital radius from the geocenter. Notice that now a relative distance variation is shown so as to describe the distance variation from the different satellite positions to the geocenter. The orbital radius is the nominal one, a circumference described by Schwarzschild metric. The orbit perturbations that are taken into account in the metric are Moon plus Sun gravitational potential. The Earth oblateness effect is not shown because it decreases when the distance from the geocenter increases. This effect is much smaller than Moon and Sun gravitational potential effects for orbital radius greater than about 5 × 10 4 Km.
From Fig. 5, we can conclude that the perturbations of a satellite orbit given by Moon plus Sun gravitational potentials increase as the satellite orbital radius increases. Fig. 6 shows the radial distance, divided by its given orbital radius, for a satellite, as a function of the proper time (for two orbital periods), at different orbital radius from the geocenter. The orbital perturbation effect that is considered in the metric is the Earth oblateness. Fig. 6 clearly shows that when the orbital radius of the satellite is increasing the Earth oblateness effect is decreasing. This effect is significant for satellites with lower orbital radius.

Relativistic positioning and user location
In this section, the U-errors are computed taking into account the difference in location when considering Schwarzschild metric and the metric we have introduced. In this way we can appreciate the difference in the positioning with the two metrics. The description of the satellite world lines taking into account the Moon, the Sun and the Earth quadrupole is more precise and is closer to the positioning it could be expected from the Galileo Satellite Constellation.
As stated before, the computations presented here are numerically calculated implementing known analytical results by Coll et al. (2010) and (2010b) concerning RPS.
Our computations are performed in a similar way as in Puchades and Sáez (2014). Let us then describe the main features of such computations. Firstly, in subsection 4.1, two codes (XT code/T Xcode ) are presented, as well as their implementation to calculate the positioning errors. Secondly, in subsection 4.2, the procedure to estimate the positioning errors is described. Thirdly, in subsection 4.3, the HEALPIx representation and the initial users distribution are explained. This kind of representation is used to display our numerical results. Fourthly, in subsection 4.4, the analysis of the numerical results is done. Finally, in subsection 4.5, the spatial configuration of user-satellites that is associated to maximum positioning errors is described.

XT code/T Xcode
The positioning of a user in space-time is done by two calculations: 1) satellite geodesics and 2) photon geodesics.
First the satellite world lines have to be calculated. Once this is done the user location is obtained in a certain space-time. It is necessary to solve the null geodesics followed by the photons in such space-time. The user inertial coordinates, x α , are defined, and the emission ones, τ A (the proper times that the user receive at the same time from the four satellites, remember that A is the number which describe each of the four satellites). In Minskowski space-time those coordinates must satisfy the following algebraic equations: where η αβ is the Minkowski diagonal matrix with η 00 = −1 and η 11 = η 22 = η 33 = 1, and the points of the satellite world lines have inertial coordinates x β A (τ A ), which must be well known functions of the proper times τ A . According to Eqs. (7), photons follow null geodesics from satellite emission to user reception. These algebraic equations may be solved by the knowledge of the satellite world lines. Also, a numerical Newton-Raphson method (Press et al. (1999)) is used. (See Puchades and Sáez (2012) and Puchades and Sáez (2014) for a more detailed description).
We can proceed in two senses. Solve Eqs. (7) for known x α and determine τ A or inversely. In the first case, Eqs. (7) may be numerically solved for the unknowns τ A by assuming that the position coordinates x α are known. Thus, the emission coordinates are obtained from the inertial ones. However, the same equations may be solved to get the unknowns x α for known emission coordinates τ A . This second case gives the inertial coordinates in terms of the emission ones (positioning); nevertheless, this second numerical solution of Eqs. (7) is not necessary since there is an analytical formula obtained in Coll et al. (2010), which gives x α in terms of τ A for photons moving in Minkowski space-time, and arbitrary satellite world lines.
In practice, a numerical code with multiple precision was designed to calculate the emission coordinates τ A (unknowns) from the inertial ones (data) by solving Eqs. (7). It is hereafter referred to as the XT-code. This code, based on the Newton-Raphson numerical method, requires the satellite world line equations; that is to say, there must be a subroutine which calculates the inertial coordinates of every satellite x α A (τ A ) for any value of τ A (see Appendices A and B). Moreover, we built up a numerical code, based on the analytical formula obtained in Coll et al. (2010), which, for given emission coordinates τ A , allows us the calculation of the user inertial coordinates x α , positioning. This code is hereafter referred to as the TX-code.

Effects of the metric on the positioning
In this subsection, the procedure to calculate the positioning errors is described.
Let us first suppose that the satellites move without uncertainties (without any perturbation effect). Their trajectories are circumferences as it corresponds to the Schwarzschild space-time (see Appendix B). These trajectories are circular orbits in the case of a spherically symmetric non rotating Earth, in the absence of external actions (nominal trajectories).
In practice, any realistic satellite world line deviates with respect to the nominal ones. If the nominal world lines are parametrized by means of their proper times, the equations of these world lines (see Appendix B) may be written as follows: x α A (τ A ) being x α A (A = 1...4) the coordinates of a given satellite A, which is a function of its proper time τ A . In the present paper, the realistic perturbed world lines are the timelike geodesics in the space-time calculated in the previous sections (see Appendix C). The satellites are orbiting at the height of the Galileo Constellation and, in the present work, only the following major effects are considered in the metric: the Earth oblateness and the gravitational effects of the Earth, the Moon and the Sun.
In the absence of deviations with respect to the nominal lines, our XT-code gives the emission coordinates τ A corresponding to any set of inertial coordinates x α and, then, from the resulting emission coordinates, the TX-code allows us to recover the initial inertial ones. The number of digits recovered measures the accuracy of our XT and TX codes. Since multiple precision is used, this accuracy is excellent.
Let us now take the above emission coordinates τ A , which are not to be varied since they are broadcasted by the satellites and received by the user without ambiguity. For these coordinates and the perturbed satellite world lines in the space-time (including the perturbations effects: Moon and Sun gravitational field and the Earth oblateness) calculated in the previous sections, the TX-code gives new inertial coordinates x α + ∆(x α ). Coordinates x α + ∆(x α ) are to be compared with the inertial coordinates x α initially assumed. The quantity is a good estimator of the positioning errors produced by the perturbations of the satellite motions. Those errors are called U-errors. It is worthwhile to emphasize that user positions x α and x α + ∆(x α ) correspond to the same emission coordinates, which are received from the satellites, but to different world lines. The nominal world lines lead to x α and the perturbed ones give x α + ∆(x α ). We may then say that the user position is x α with an error whose amplitude is given by the estimator ∆ d . See Puchades and Sáez (2014) for a more detailed explanation of the definition of these U-errors. The improvement presented here is the use of a most accurate description of satellite perturbations using a metric which better accounts of a more accurate trajectory of the satellites.

HEALPIx representations and initial users distribution
In a similar way as it was performed in Puchades and Sáez (2014), we are going to present the results of the computation of the ∆ d , the U-errors, of different users equidistributed in sphere surfaces with different radius centred in the geocenter. In order to represent such values an appropriate method is needed. This method allows us to represent some quantities in 3D, t = constant, space-time sections. An appropriate pixelization and colour bars are needed.
The HEALPIx (Hierarchical Equal Area Isolatitude Pixelisation of the Sphere) package (Górski et al. (1999)) is a very useful tool for this representation and it is used here. This package was initially designed to represent the Cosmic Microwave Background temperature distribution in the sky. As it writes the values of a scalar quantity in a sphere surface, it is also very convenient for our kind of pictures. This method displays any quantity as a function of the direction (pixel). The sphere surface is divided in 12 × N 2 side pixels and the free parameter N side takes even natural numbers. In order to better compare with the pixelization considered in Puchades and Sáez (2014) N side = 16 is taken. Initial users are fixed on a surface of a sphere of radius R centred in the geocenter. One user per each HEALPIx direction. 3072 initial users (pixels) equally distributed are considered. The angular area of any one of them is ≈ 13.43 square degrees. Such angular area is close to sixty four times the mean angular area of the full moon. However their shape in the mollweide representation used here is not the same for all pixels. They are more elongated in the polar zone. For each initial user, we carry out the procedure that we have commented for the computations. The pixelized sphere is shown by using the mollweide projection, in which, the frontal hemisphere is projected on the central part of the figure, and the opposite hemisphere is represented in the lateral parts.
As stated in the latter paragraph, HEALPIx mollweide maps are taken to represent the scalar quantity required, in our case the U-errors, ∆ d . According to the colour bar displayed, any pixel shows a colour which states the ∆ d value in the map for the direction associated to such pixel. ∆ d values for different cases are shown in Figs. 7, 8 and 9, considering Galileo satellites 2, 5, 20 and 23, in a similar way as in Puchades and Sáez (2014). HEALPIx mollweide maps of Figs. 7 and 8 are ∆ d values, in Km, on spheres with Earth radius 6378 = R ⊕ at different user times t. Otherwise, Fig. 9 shows ∆ d values, in Km, on spheres with different radius at user time t = 19 h. From top to bottom, the radius of the spheres in kilometres are 1.5 × 10 4 , 3 × 10 4 and 5 × 10 4 . The first and the third values are two of the values taken in Puchades and Sáez (2014). This fact allows us to compare the results obtained with the two different procedures to compute the U-errors.

Analysis of results (Figs. 7 and 8)
In this subsection, the analysis of the numerical results for different initial users distributions is done. As stated at the end of the previous section, sphere surfaces at the Earth radius 6378 = R ⊕ are taken. We want to know which is the U-error, ∆ d , produced on the positioning taking into account the deviation caused on nominal orbits (Schwarzschild ones) when the effects of Moon, Sun and Earth quadrupole are considered in the metric. In Puchades and Sáez (2014) it was pointed out that for 6378 = R ⊕ sphere radius the positioning errors where of the same order as the satellite deviations from the nominal orbits. In the cited paper, the deviation was statistically generated. Now more realistic perturbations are considered. Satellites with different relative positions are considered at different times.
The results commented in the last paragraph are obtained from HEALPIx maps shown in Figs. 7 and 8. Let us make a list of our main conclusions: • The greatest ∆ d values (red pixels) correspond to having the maximum radial distance deviation of the satellite for the case of the four chosen satellites (see Fig. 2, violet line). When the satellite trajectory deviation with respect to the Schwarzschild trajectory (considering circular orbit as nominal trajectory) increases, then the value ∆ d increases. For the positioning, the four satellites are needed. Therefore, this deviation for the four chosen satellites has influence on the positioning errors. So, the value ∆ d depends directly on the satellite-Earth-Moon-Sun spatial configuration, and this for the four satellites.
• The ∆ d values are of the same order of magnitude of the perturbation applied on the satellite orbit in most of the pixels.
That is, about a few Kilometres is this order of magnitude, just the same as the satellite deviations when we take into account Earth oblateness plus Sun plus Moon gravitational field. These results are in agreement with the respective results obtained in Puchades and Sáez (2014).
• Almost the same shape on the maps (but with different maximum and minimum values) is repeated when we consider the next satellite orbital period. This periodic effect was to be expected because the Sun-Moon-Earth-satellite spatial configuration practically does not change after one satellite orbital period (14.2 h). Therefore, Sun, Moon and Earth influence on the satellite orbit is almost the same than in the previous satellite orbital period. This periodic effect is shown in Figs. 7 and 8, corresponding to the cases where the receiver time is 11h and 25h, respectively.
• The Earth oblateness (Earth quadrupole) effect is shown in Fig. 8 (receiver time = 21 h), as it can be seen in Eq. (5), the sinusoidal shape is obtained.

Jacobian and user-satellites spatial configuration (Fig. 9)
In this subsection we present a case with a user-satellites spatial configuration that is associated to maximum positioning errors. The purpose of it is to enhance the relationship of the U-errors with the Jacobian, J, values as it was stated in the Puchades and Sáez (2014) work. For our computation of the errors based on a more realistic case, the relationship between the Jacobian values near zero and the great values of the U-errors still stands, as it must be. For the sake of comparison the fig 2 top case of Puchades and Sáez (2014) has been considered. Let us recall which is this Galileo satellites configuration: satellites 2, 5, 20 and 23 at user time t = 19 h. Fig. 9 shows HEALPIx mollweide maps of ∆ d values, in Km, on spheres with different radius for this satellite distribution. Two radius spheres in kilometres as in figures 5 and 6 of Puchades and Sáez (2014) have been considered: 1.5 × 10 4 (top) and 5 × 10 4 (bottom). Also 3 × 10 4 (middle) is taken here. Top subfigure can be compared with the top one of fig 6 in Puchades and Sáez (2014), the satellites configuration and user time are the same ones. The difference between both figures is the way of computing the ∆ d errors. The distribution of those U-errors is very similar as it can be seen comparing both subfigures. This is logical as the same satellites configuration and user time have been taken. When we use our metric, lower errors are obtained. The range in Puchades and Sáez (2014) varies from 6.84 to 24.0 km, while in the present paper the range of such values varies from 0.2 to 10.1 km. The lower values obtained in the present work should be due to our better approximation: taking a realistic metric better describes the U-errors in RPS.
As in the mentioned paper, some cutoffs of ∆ d are considered. In Fig. 9, grey coloured pixels are characterized by the condition ∆ d > 50 Km for the middle subfigure and ∆ d > 100 Km for the bottom subfigure. The same as in our previous work, as the radius of the sphere increases the ∆ d values increase (more grey pixels). In the above cited paper it was explained the relation between the volume of the tetrahedron formed by the tips of the four user-satellites unit vectors, V T , and the Jacobian, J. Specifically, this volume is a sixth of the |J| value: V T = |J|/6 (see also Langley (1999)).
We have chosen a user whose ∆ d value is one of the greatest of Fig. 9 (R = 3 × 10 4 Km). Fig. 10 shows the locations of the four satellites in the celestial sphere of this user, whose user-satellites spatial configuration corresponds to a case of J ≃ 0 (J = −6.7 × 10 −2 ). Let us recall that J = |∂τ A /∂x α |, of the transformation giving the emission coordinates in terms of the inertial ones. Note that the considered Galileo satellites are 2, 5, 20 and 23 at user time t = 19 h for a user on a surface of a sphere of radius R = 3 × 10 4 Km. Remember that the volume V T is proportional to the Jacobian value. In Fig. 10, the satellites are represented by four black asteriks, the user by a blue cross (the origin of the Cartesian system) and the centres of the circles by red, blue, green and purple asterisks. Four different circles, that pass through three points, with each possible combination of three satellites (colour points) are also drawn. As it was explained in Puchades and Sáez (2014), such satellite configuration corresponds to a tetrahedon near of null volume value and therefore to a J close to zero. This fact accounts for such a great U-error value.
In Puchades and Sáez (2014) it was also studied the relation of the user distance to Earth and the ∆ d errors. Let us make a comment here taken into account the figures presented in this work. From HEALPIx maps shown in Fig.9, we also see that as the radius of the spherical surfaces increases the ∆ d values also increase (see how the minimum value of ∆ d increases with R in these maps). This is because if the user is very far from the satellites, they are all in a small solid angle and the tetrahedron volume V T is also expected to be small. Moreover, on these maps, there are regions with J ≃ 0 (grey pixels with high ∆ d value) when the height is greater than approx. 2 × 10 4 Km. In Puchades and Sáez (2014), it was also shown that not only regions near J null values give great ∆ d error values. There are also regions where ∆ d is small. One of the conclusions in the mentioned paper was that to solve this problem of big positioning errors at great distances, it should be interesting in the future, to choose the best combination of four satellites (from GPS and Galileo constellations), whose user-satellites spatial configuration gives smaller ∆ d errors. This fact seems possible if the satellites configurations taken into account are changed accordingly. Some proposals of how to determine such configurations can be seen in Puchades and Sáez (2014) and also apply to the results obtained in this paper. Another way to do that could be to locate the satellites in other orbits in the solar system, in a way that the user-satellites spatial configuration corresponds to V T values not close to zero.

Conclusions and perspectives
The main purpose of this paper is to study the U-errors obtained as the difference in RPS by using Schwarzschild and a more accurate metric to describe the Galileo satellite world lines. This represents an advance with respect to the research made by Puchades and Sáez (2014). In Puchades and Sáez (2014) the satellite perturbations were statistically computed. Here such perturbations are computed differently by taking into account the effects of Moon, Sun and Earth quadrupole in the metric. These are the greatest effects perturbing satellite world lines at the height of Galileo or GPS constellation. In order to better understand the difference in the positioning, the U-errors, a deep study of the change in the Galileo satellites trajectories has been first developed, see sections 2 and 3. This study takes into account the contribution on the satellite world lines of such three terms by (i) taking them separately, (ii) combining two of them and (iii) considering all three together.
A Runge-Kutta algorithm is used to solve the geodesic equations. High accuracy is achieved (10 −18 ). Also multiple precision (40 significant digits for each number) is used. Adequate initial conditions have been found to solve the ODE. Precise satellite geodesics are required to be able to precisely compute the user location in space-time. Moreover, this precision is needed to incorporate other small contributions to the satellite orbits perturbations. In this work, only the Galileo Satellite Constellation is considered. As it is known, one of the greater contributions at the height of Galileo satellites is the Earth oblateness (about a few kilometres). The Moon and Sun gravitational effects are also important (about a few hundred Kilometres). Moreover, due to the greater relative Moon motion, after one orbital period, the satellite position is shifted when the Moon is considered. Our results are in agreement with other results known in the literature. Therefore, our method is working properly because the order of magnitude computed in preceding works with the GNSS approximation is obtained (see for instance Teunissen and Montenbruck (2015), Chapter 3, Montenbruck and Gill (2005) and Roh (2018)). The RPS methods are more exact than GNSS classical procedures, even when the last ones incorporate relativity corrections. From the beginning, we numerically solve the satellite geodesic equations (ODE) in our given space-time.
Then we have simulated satellite world lines at different distances from Earth (see section 3) and studied the influence of the perturbations considered in our metric. The Earth oblateness orbit perturbing contribution is of the same order of magnitude than that of Moon and Sun gravitational effect at 5 × 10 4 Km. As the satellite orbital height increases, the Earth oblateness effect decreases and Moon and Sun effects increase. At 10 5 Km, the Earth oblateness effect is smaller than the Moon and the Sun contributions. At these distances, the satellite position is also shifted when the Moon presence is considered, as it happens at the Galileo satellites distances. This is an improvement with respect to our previous works because now it is possible to appreciate the separate contribution of each perturbation in the RPS.
After studying the satellite trajectories, we carry out RPS with our metric. A similar procedure as in Puchades and Sáez (2014) has been performed. But now, a metric taking into account the greater physical perturbations at Galileo Satellite Constellation is considered to compute ∆ d . Notice that the same formula, Coll et al. (2010), has been used to compute the proper times that the user receives from the satellites. Also HEALPIx maps are considered to describe the U-errors as in Puchades and Sáez (2014).
As it was concluded in the cited paper, the positioning errors values, ∆ d , are almost of the same order of magnitude as those of the perturbed satellite orbits (orbital perturbation effect). Here the highest ∆ d values (red pixels in subsection 4.4) correspond to having the maximum radial distance deviations of the satellite for the case of the four chosen satellites in Fig.  2 (violet line). So, the value ∆ d depends directly on the satellite-Earth-Moon-Sun relative spatial configuration and it does so for each of the four satellites. Almost the same HEALPIx maps are recovered after a Galileo satellite orbital period. This is because the relative spatial configuration among satellite-Moon-Sun-Earth does not nearly change after 14.2 h (periodic effect), as the Moon and Sun hardly move after a Galileo orbital period. Here the ∆ d values are smaller than the ones obtained with the statistical procedure used in Puchades and Sáez (2014). This fact is due to a more realistic representation of the satellite orbital perturbations by the use of metrics instead of statistical deviations (see subsection 4.5). Therefore a more accurate computation of the U-errors is performed and so a more precise calculation of the user's positioning can be achieved.
As in Puchades and Sáez (2014) when the radii of the spherical surfaces increase, the ∆ d values increase. This is because if the user is very far from the satellites, these are all in a small solid angle and the tetrahedron volume V T is expected to be small. Moreover, on these maps, there are regions with J ≃ 0 when the radius is greater than about 20000Km. One solution when J ≃ 0, is to choose another combination of four satellites from Galileo or GPS constellation, whose user-satellites configuration is associated to values of J not close to zero. An alternative way to do that could be to locate the satellites in other orbits in the solar system, in a way that the user-satellites spatial configuration corresponds to V T values not close to zero.
We are currently working on an improvement in our numerical procedure. The Newton-Raphson numerical method (see XTcode in section 4.1) could be avoided, in such a way that analytical functions of the proper time are not used, since they imply the use of Schwarzschild world lines for the satellites. For example, to use the secant method with the satellite world lines (taking into account the Sun and Moon presence and the Earth quadrupole). We think this will improve the numerical code and results.
Once this is done, another interesting thing to do, and very useful, should be to create HEALPIx maps as in subsection 4.4 (see Figs. 7 and 8), but calculating the positioning error ∆ d on the geoid, instead of on the spherical surface of radius R ⊕ . It would allow us to calculate the accuracy of Galileo satellites, with our implementation, and compare with the data from Galileo Constellation, and other constellations. These results should also be interesting for geodesic treatment.
There are other perturbations a part from those considered in this paper that also contribute to the computation of the satellite world lines. The order of magnitude of such contributions depends on the satellite's height as it can be seen, for instance, in Montenbruck and Gill (2005) (see, in particular, fig 3.1 at page 55). Therefore, this variation of the trajectories of the four satellites considered should contribute to the change in the calculations of RPS. We are studying such contributions and will present the results obtained elsewhere.
Also the use of our method in space navigation is being planned. The Barycentric Celestial Reference System could be used as reference system to locate the emitters (four satellites) in the solar system, in such a way that the configurations of the user-satellites associated to J ≃ 0 are avoided. For example, in the vicinity of the Moon, two emitters fixed on the Moon surface (North and South poles) and two emitters from Galileo satellites. In such a case, the V T is not close to zero, in most of the cases. The positioning of a spacecraft that navigates in the solar system could be determined considering emitters in other appropriate locations to be studied. Teunissen, P., Montenbruck, O., 2015. Global navigation satellite systems, Springer Handbook, 1st ed. Cham, Switzerland.

A Effects on the Galileo satellite world lines
So as to better understand the perturbations on the satellite world lines, at the beginning of our work we made some computations of the order of magnitude of the effects of the Earth, Sun and Moon on Galileo satellites. These are the major perturbations on the satellite orbits at the satellite height (see Roh (2018)). For this reason it is interesting to study the magnitude of their gravitational field and potential at approx. 3 × 10 4 Km from the geocentre (the radius of a Galileo orbit). Also the relation between them (Earth-Sun and Earth-Moon system) is relevant. In a near future the order of magnitude of other smaller effects will be studied. Notice that here we do not present the metric itself, we only compare the relation between the order of magnitude of Newtonian effects on satellite orbits, for the sake of better understanding the impact of the effects considered here. We assume that the geocentric position vector for Sun, Moon and satellite is x ⊙ = (−1, 0, 0) AU , x = (−384402.0, 0, 0) Km and x = (29655.3, 0, 0) Km, respectively. In all the paper ⊙, ⊕, stands for Sun, Earth and Moon respectively. As it has been said the origin of the reference system is the geocentre (GRS). Afterwards, when we consider the Earth movement and acceleration, the GRS should be co-moving (the origin of the GRS is moving with the Earth centre).
The quantities computed here are the order of magnitude of Newtonian potentials and accelerations on the Galileo satellites positions.
Two cases are computed: (a) Stationary Earth.
(a) Earth with accelerated motion and the reference system co-moving with the Earth.
For the (a) case, we calculate: • The Earth potential divided by the Sun potential: φ⊕ φ⊙ = 0.0150. The result is that the Sun potential is 67 times greater than the Earth's potential.
• The intensity of the Earth gravitational field divided by the intensity of the Sun's gravitational field is: a⊕ a⊙ = 74.6525. Then, the force of the Sun on the satellite is 75 times smaller than the force of the Earth on it.
Although the Sun potential has a great effect, its variations, a ⊙ , corresponding to the forces which represent the gravitational field, can be neglected in a ⊕ and the quantity φ ⊙ could be taken as an additive constant of the φ ⊕ .
• Moon potential divided by the Earth's potential: φ φ⊕ = 0.0010. This effect (φ ) is either negligible or could be taken as a perturbation.
• The intensity of the Moon gravitational field divided by the intensity of the Earth's: a a⊕ = 8.8139 × 10 −5 . This quantity (a ) could be neglected (as we do) or considered as a small perturbation.
The Newtonian concept of force cannot be applied to GR. The Christoffel symbols take the "role" of gravitational forces. In the set of ODE, Eqs. (30)-(33), as it can be seen in Appendix C, the Sun, Earth and Moon "gravitational potentials" are placed in the denominator. And the term 1 + ∂φ ⊕ + ∂φ ⊙ + ∂φ is in the numerator. We next present a summary of the results obtained for case (a): Results of case (a): However, the Earth is accelerated [(b) case] under the action of the celestial body (Sun or Moon) and such action will be included now.
M B is the mass of the celestial body B (Sun or Moon), x B its geocentric position and x is the geocentric position vector of the satellite.
The term a B1 of Eq. (9) is: This term a B1 stands for the perturbing acceleration of B exerted on the satellite. And the term a B2 of Eq. (9) is: This term a B2 is the perturbing acceleration (on the satellite) caused by the action of B on the Earth. This is the inertial contribution (notice the minus sign) caused by the accelerated GRS. Accelerations (a B1 and a B2 ) on the orbit of a Galileo satellite due to the Sun and the Moon are computed and their values can be seen in Table 1 Table 1 shows the first Cartesian component of a B . The numerical results of (b) appear in the fourth column of Table 1. Nevertheless, when we consider a stationary Earth [(a) case], the term a B2 is zero. The results in case (a) appear in the second column of Table 1.
We can summarize the results obtained and presented in Table 1 in this way: i) the perturbing acceleration from the Sun is slightly smaller than that of the Moon. ii) Meanwhile, the terms a B1 and a B2 from the Sun are 100 times greater than those from the Moon. iii) As a first conclusion, we observe that the intensity of the Sun gravitational field is smaller than such from the Moon at the GRS that considers the Earth accelerated by the gravitational force of the Sun and the Moon.
On the one hand, a stationary Earth (in a GRS), gives an acceleration from the Sun 100 times greater than that from the Moon. On the other hand, we have obtained an acceleration from the Earth central force ( The effect of the oblateness of the Earth is now studied. Its Cartesian components (a J2 1 , a J2 2 , a J2 3 ) with regard to its acceleration, (Seeber 2003, Sharma et al. 2019, are given by where A(i) is equal to 1 for i = 1, 2 and A(i) is equal to 3 for i = 3. The perturbing acceleration from the Earth oblateness (for the first Cartesian component), a J2 1 , has been computed and its value is −6.793 × 10 −5 (m/s 2 ) for θ = 0, as it is shown in the third row of Table 1. Let us recall that the satellite position vector is x = (29655.3, 0, 0) Km.
As it can be seen in the fourth column of Table 1, the Earth oblateness produces an acceleration on the orbits of the Galileo satellite one order of magnitude greater, in absolute value (6.793 × 10 −5 ) , than those produced by the Moon potential (4.583 × 10 −6 ) or the Sun potential (2.350 × 10 −6 ).
In this Appendix, we have stated some important relations. They have helped us to better interpret the results obtained in our work.
B A short description of the circular satellite world lines As it is well known, in the unperturbed Minkowskian space-time, the spatial location of a Galileo satellite A, which moves along a circumference, requires three angles. Two of them, Θ and ψ, characterize one of the three orbital planes. These angles are constant. The third angle, α A , localizes the satellite on its trajectory. It depends on time. All this was taken into account in Sáez (2012, 2014) to find the world line equations of the satellite A to first order in the small dimensionless parameter GM ⊕ /R, whose maximum value is GM ⊕ /R ⊕ ≃ 6.94 × 10 −10 . These equations are as follows: where the factor γ and the angle α A are given by the relations (Ashby 2003, Pascual-Sánchez 2007 and  (2000), Soffel et al. (2003)) recommends a resolution for the metric tensor in the Geocenter Celestial Reference System GCRS up to order o(1/c 2 ) (neglecting terms smaller than o(1/c 2 ), for the A o(1/c 3 ) body vector potential), such as: where (t ≡ geocentric coordinate time, x) are the pseudo-Cartesian isotropic GCRS coordinates, w 0 = G A M A /r A summation over all solar system bodies, r A = x − x A , x A is the position vector of the centre of mass of A body, r A = |r A |, and where w L has the expansion in terms of multipole moments that requires each body. In this paper only the masses M A of three bodies are considered: Sun, Earth and Moon (A = ⊙, ⊕, ). Although there are neglected terms with order greater than GMA rA , they will be considered in the future. w L has been considered, being 2w L = φ J2 and φ J2 the Earth quadrupole potential. The quantity w 0 as w 0 = G M⊕ r⊕ + M⊙ r⊙ + M r is also taken. φ is defined as 2 (w 0 (t, x) + w L (t, x)).
A numerical integration of timelike geodesic equations is performed to compute the satellite trajectories x α (τ ) in the GCRS, where u ν is the four-velocity (u ν = dx ν dτ ) and τ the proper time. The covariant and contravariant components of the metric in pseudo-Cartesian isotropic GCRS coordinates have this form: and the Christoffel symbols are calculated from this definition: with: being ζ ,α ≡ ∂ζ ∂x α , and ζ a function of x α .
The potential φ = φ(x i ) is considered as stationary. The non-null Christoffel symbols are then obtained: The Christoffel symbols Γ k ij are as follows: (no summation over repeated indices).
The following expressions are the geodesic equations for these Christoffel symbols: The system has the following constraint: This constraint lets us verify the proper functioning of the code. In each step of the numerical integration the code checks that Eq. (34) holds. A Runge-Kutta method is applied to solve such equations. In order to completely describe the motion of the system considered here, the motion equations of the celestial bodies taken into account should also be added to the system of Eqs. (30)-(33). The Newtonian equations are sufficient. The Newtonian Earth motion equations in the GCRS, for the system Earth-Moon-Sun, are given by: similar equations for the Sun stand: and also for the Moon: The greatest φ "potential" contribution is the Earth's contribution, which has the following Schwarzschild metric form (see Appendix A): where r ⊕ = x − x ⊕ . Additionally, the satellite geodesics in Schwarzschild perturbed space-time need also be computed. The gravitational potential is therefore: where φ pert are the additional perturbing potentials produced by the Sun, the Moon and the Earth quadrupole. The quantities φ ,i are the first spatial derivatives ∂φ/∂x i , where x i is a function of the proper time x i = x i (τ ), so that: These are the spatial derivatives of the Earth's gravitational potential (central force): It should be noted that in this case the reference system is co-moving with the Earth (GRS). In order to calculate the satellite world lines the time-like geodesic equations from this metric are numerically solved. To solve the ODE system presented here a Runge-Kutta method with adaptive step has been implemented. Our numerical algorithm has high accuracy (10 −18 ) and multiple precision (40 significant digits for each number). This precision will also allow us to consider other PN effects as our work advances. Our algorithm is being prepared to provide the precision and accuracy needed for these future additions. Notice once more that the effects computed in the present paper consider the corresponding terms of Moon or Sun or Earth oblateness or a combination of them.  Earth oblateness (green) and Moon plus Sun effect plus Earth oblateness (violet) are here taken into account. Fig. 3: Radial distance of a satellite versus the proper time (for two orbital periods) at 5 × 10 4 Km from the geocenter. The orbital perturbation effects that are considered in the metric are Earth oblateness (green), Moon gravitational potential (blue), Sun gravitational potential (red), Moon plus Sun gravitational potential (magenta) and, Earth oblateness plus Moon plus Sun gravitational potential (violet). Fig. 4: Radial distance of a satellite versus the proper time (for two orbital periods) at different orbital radius from the geocenter. From top to bottom, the considered orbital radius is: 10 5 and 1.5 × 10 5 in Km. The orbital perturbation effects that are considered in the metric are Moon gravitational potential (blue), Sun gravitational potential (red), Moon plus Sun gravitational potential (magenta) and, Earth oblateness plus Moon plus Sun gravitational potential (violet).  Fig. 6: Radial distance R, divided by its nominal orbital radius, R 0 , that is shown in the legend, for a satellite, as a function of the proper time (for two orbital periods), at different orbital radius from the geocenter. The orbital perturbation effect that is considered in the metric is the Earth oblateness (green). h GAL means the height of a Galileo satelite.    10: Locations of the satellites in the celestial sphere of a user, whose user-satellites configuration corresponds to a case of J ≃ 0. The considered satellites are 2, 5, 20 and 23 at user time t = 19 h for a user on a surface of a sphere of radius R = 3 × 10 4 Km. The satellites considered are represented by four black asterisks, the user by a blue cross (the origin of the Cartesian system) and the centres of the circles by red, blue, green and purple dots. Four different circles, that pass through three points, represent each possible combination of the three satellites (points).