A Simplified–Modified Algorithm for Glonass Broadcast Orbits Computation

The Global Navigation Satellite System (GNSS) is the worldwide set of satellite navigation constellations. Nowadays, there are four GNSSs accessible for personal purposes and free of charge (Oliveira et al., 2018). The American Global Positioning System (GPS) and the Russian Global Navigation Satellite System (GLONASS) were the first operational GNSS systems, followed by Galileo, the European system, named after the astronomer Galileo Galilei (1564–1642), and BeiDou, the Chinese system, named for the Chinese pronunciation of the Big Dipper constellation (Petrovski, 2014).

After many years of decline, restoration and modernization started in 2001; Russia’s system returned to full operation since the end of 2011 (Petrovski, 2014; Polischuk et al., 2002; Sarkar and Bose, 2017) and currently, more than 24 satellites are operational. The GLONASS system, or more precisely, GLObal’naya NAvigatsionnaya Sputnikovaya Sistema in Russian (Habrich, 1999; Shi and Wei, 2019), is like the other GNSS systems, a radio-navigation system very similar to GPS in general concept (Levine, 2002; Pace et al., 1995) that provides positioning and timing information for the users. GLONASS is operated by the Russian Space Forces (Habrich, 1999).

Coordinate and time reference frames are both key parameters of a satellite navigation system (Chen et al., 2013). In terms of coordinate reference, the Russian GLONASS system uses a reference frame called Parametry Zemli 1990 or PZ-90 (Levine, 2002) based on the PZ-90 ellipsoid; the fundamental geodetic constants and other significant parameters of PZ-90 are given in ICD-GLONASS (2016). As the time reference frame, GLONASS is based on GLONASS Time (GLONASST) synchronized to the Coordinated Universal Time UTC (SU) of Russia (Gurtner, 2015; Shi and Wei, 2019).

According to the GLONASS constellation status given in the Information Analytical Centre (IAC), official website (https://glonass-iac.ru/en/sostavOG/) of the Russian Federal Space Agency (Sarkar and Bose, 2017), in September 2024, the GLONASS space segment contains a constellation of 26 satellites in orbit, with 24 operational and two in the flight tests phase (GLONASS numbers 1 and 3). These satellites are equally distributed in three orbital planes and separated by 120° along the equator (Pace et al., 1995). The revolution period of GLONASS satellites is 11 h, 15 min, and 44 s, repeating the geometry every eight sidereal days (Petrovski, 2014; Sanz et al., 2013), and the orbital velocity is about 3.95 km/s (Levine, 2002). A review of GLONASS constellations is given in Sarkar and Bose (2017).

The ephemerides transmitted by GLONASS satellites describe the position, velocity, and luni-solar (L–S) perturbation of each satellite. This information is used to determine the coordinates of GLONASS users (receivers) by the trilateration or triangulation principle (Oliveira et al., 2018) using the known coordinates of GLONASS satellites, which can be computed at any time by the Runge–Kutta (R–K) integration method, which is the subject of this paper.

Several works have used the R–K method in the computation of GLONASS orbits; a literature review is given in Maciuk (2016). In Lin et al. (2009), the authors used the fourth order, where the L–S acceleration has been taken as constant, and the fourth and fifth orders, where the L–S acceleration has been taken as a variable. In Maciuk (2016), the author used several orders (fourth, fifth, Fehlberg fourth, and Fehlberg fifth) of the R–K method. In Góral and Skorupa (2015), the authors present the analytical method to interpolate and predict the position of GLONASS satellites where the total lunar–solar accelerations were calculated according to the formulae given in the ICD-GLONASS document. In Medjahed et al. (2021), the L–S acceleration was included in the order of 04 of the R–K method.

The GLONASS Interface Control Document (ICD) gives the precise and simplified algorithms used in the computation of the broadcast orbits (position and velocity) of GLONASS satellites using the R–K integration method. These algorithms are defined according to the integration interval and L–S taken into account during the integration period of the Ordinary Differential Equation (ODE) describing the motion of GLONASS satellites as follows (ICD-GLONASS, 2016):

Precise algorithm (Algorithm J-1, page 45): The L–S acceleration terms are computed according to the formulas given in this algorithm.

In both algorithms, the differential equation to be resolved by the R–K method is composed of six sub-equations with six initial conditions.

In this study, the simplified algorithm is modified and called the simplified–modified algorithm, where the three L–S accelerations have been added as variable terms into the equation of GLONASS satellites' motion. The derivatives of these accelerations required for the integration by the R–K method have been determined from linear functions recalculated every 30 min, the period of GLONASS broadcast ephemerides validity. The system of equations becomes composed of nine sub-equations resolved by the MATLAB ode45 solver to get the position, velocity, and L–S acceleration of GLONASS satellites at any desired time from the broadcast ephemerides used as initial conditions.

This paper is organized as follows: in Section 1, we have given a general presentation of the GLONASS system and the different algorithms used in GLONASS orbit computation. In Section 2, the ODE of GLONASS satellite motion and the linear functions computation along each axis are given; in Section 3, we present the ode45 solver implementation in MATLAB. In Section 4, the data source used is presented first, followed by a comparison between the simplified algorithm, that is, the L–S acceleration taken as a constant, and the simplified modified algorithm, that is, the L–S acceleration taken as variable terms. Finally, a summary of the conclusion is given.

The motion of the GLONASS satellites orbiting the Earth, including the L–S acceleration, is given by the ODE of second order transformed into first order given in Lin and Guo (2009), Maciuk (2016), Medjahed et al. (2021), and Sanz et al. (2013) by the following equation: 1{ Y(˙1)=Y(4)Y(˙2)=Y(5)Y(˙3)=Y(6)Y(˙4)=−GMr3Y(1)+(32)C20GMa2r5(1−5Y(3)2r2)Y(1)+w2X+2we2Vy+Y(7)Y(˙5)=−GMr3Y(2)+(32)C20GMa2r5(1−5Y(3)2r2)Y(2)+w2Y2−2we2Vx+Y(8)Y(˙6)=−GMr3Y(3)+(32)C20GMa2r5(3−5Y(4)2r2)Y(3)+Y(9)Y(˙7)=axY(˙8)=ayY(˙9)=az

The coefficients of Equation 1 are given in the ICD-GLONASS (2016) document as follows: GM is the mass of Earth; r is the orbital radius; C₂₀ is the gravitational constant; a is the semi-major axis of the PZ-90 ellipsoid, and w is the rotation rate of Earth. The value of each coefficient is given in Appendix 1.

The coefficients y are the ephemerides of satellites given in the broadcast file and used as initial conditions where Y(1), Y(2), and Y(3) are the X, Y, and Z positions; Y(4), Y(5), and Y(6) are the X, Y, and Z velocities (km/s), and Y(7), Y(8), and Y(9) are the X, Y, and Z L–S accelerations of GLONASS satellites (km/s²).

a_x, a_y and a_z are the derivatives of the linear functions describing the variation of the L–S acceleration between two consecutive epochs. This variation along the X-, Y-, and Z-axes is given by: 2{ Y(7)=ax·t+bx=>Y(˙7)=axY(8)=ay·t+by=>Y(˙8)=ayY(9)=az·t+bz=>Y(˙9)=az

Where (t) is the variable time t = [t_b t_b + 30 min] and (t_b) is the time of ephemerides, where t_b = hour + 15 min or t_b = hour + 45 min.

In the document Cook (1962), the author began the introduction with the following sentence: ‘‘The effects of the gravitational attractions of the Sun and Moon on the orbits of artificial Earth satellites are now important for two reasons….” The L–S acceleration describes the combined gravitational effects of moon and the sun on the motion of an artificial satellite as given in Figure 1. A part of the research published in the 50s of the last century on the effects of moon and sun on the artificial satellites is presented in Cook 1962.

Figure 1.

Moon and sun perturbations

Several papers have examined the impact of the sun and moon perturbations on the artificial satellites. In Solórzano and Almeida Prado (2013), an overview of the most recent publications is included.

The ODE of satellite motion can be solved analytically or numerically (Maciuk, 2016). The numerical methods find approximate solutions to the complicated model equations (Banu et al. 2021). There are numerous numerical methods used to solve differential equations (Ardourel and Jebeile, 2016). The Euler method is one of the simplest and oldest numerical methods proposed by Leonhard Euler in 1768 (Banu et al., 2021). The R–K methods are also used to solve ODEs; these methods are particularly easy and may be applied to a wide range of different problems (Montenbruck and Gill, 2012). The Bulirsch–Stoer method is a highly accurate extrapolation method that combines solutions from different step sizes to reduce the error; it is useful when a highly accurate solution is required with minimal computational effort (Press et al., 2007). The leapfrog method is commonly applied in physics simulations due to its ability to conserve energy in long-term simulations. The Dormand–Prince method adjusts the step size dynamically to maintain a desired error tolerance (Mathworks, 2024).

Each method has its advantages, disadvantages, and area of application. In this paper, we use the Dormand–Prince R–K method adapted in the MATLAB ode45 solver to resolve the ODE of GLONASS satellite motion given by Equation 1; ode45 is one of the most popular solvers used to solve differential equations (Waleed, 2013) given in the documentation by the following equations (Butcher, 2016; Dormand and Prince, 1980; Press et al., 2007): 3Yi+1=Yi+(35384K1+0.K2+5001113K3+125192K4−21876784K5+1184K6)

To get the next value Y_i+1 from the previous value Y_i (i.e., position, velocity, and L–S acceleration), there are six k-coefficients for each component given by: 4{ K1=h.f(Xn,Yn)K2=h.f(Xn+h5,Yn+K15)K3=h.f(Xn+3h10,Yn+3K140+9K240)K4=h.f(Xn+5h4,Yn+44K145−56K215+32K39)K5=h.f(Xn+8h9,Yn+19732K16561−25360K22187+64448K36561−212K4729)K6=h.f(Xn+h,Yn+9017K13168−355K233+46732K35247+49K4176−5103K518656)

The following syntax implements the ode45 solver for the numerical resolution of Equation 1 (Chapra, 2008; Mathworks, 2024): [Output] = Solver (Input) where:

[Output]=[t, Y] where:

Broadcast ephemerides describe the orbital motion through a set of parameters over a limited validity range of less than a few hours (Montenbruck et al., 2020). For the GLONASS system, the validity of these parameters is 30 min given in the PZ-90 reference frame (Góral and Skorupa, 2015; Lin and Guo, 2009).

The files containing this ephemeris, calculated by the ground control segment, include a header giving general information, then a list of parameters, such as positions, velocities, and L–S accelerations of each satellite. The description of these files can be found in Gurtner (2015), Maciuk (2016), and Medjahed et al. (2021).

In this study, we computed the broadcast orbits (positions, velocities, and L–S acceleration) of 24 GLONASS satellites between March 1 and 21, 2024 (21 days).These data are available to the scientific community through the global data center Crustal Dynamics Data Information System (CDDIS) located at the National Aeronautics and Space Administration or NASA’s Goddard Space Flight Center in Greenbelt, MD, USA (https://cddis.nasa.gov/About/Background.html).

Figure 2 shows the orbit of GLONASS satellite number 1 given in the data file (black dots) and the orbit computed by the ode45 solver (continuous magenta lines).

Figure 2.

Orbits of GLONASS satellite number 01 (March 8, 2024)

The ode45 solver was used to resolve Equation 1 using the simplified and the simplified–modified algorithms.

Table 1 summarizes the statistics of differences (min, max, mean, root mean square [RMS]) in position and velocities on the X, Y, and Z axes of 24,168 epochs of computation (24 satellites, each computed at time H: 15 min and H: 45 min during 21 days). The results of these differences are presented as follows:

Simplified algorithm (Δ₁): The difference between the broadcast ephemerides provided in the data file noted Y₀ and the integrated solution using the simplified algorithm noted Y_ode45(1), where the L–S accelerations are taken as constant. Δ₁ = |Y₀ – Y_ode45(1)|

Figures 3 and 4 show the differences between the results of the simplified–modified algorithm and the broadcast ephemerides data.

Figure 3.

Difference (Δ₂) in X, Y, and Z (positions)

Figure 4.

Difference (Δ₂) in Vx, Vy, and Vz (velocities)

Figure 4 shows the differences in position, velocity, and L–S acceleration between the algorithms used in this study (i.e., simplified algorithms and simplified–modified).

The maximum difference values: Δ= |Δ₁ - Δ₂| on the X, Y, and Z axes are given in Table 2.

The 3D maximum difference in position is 1.76 m, in velocity is 0.0027 m/s, and in the L–S acceleration is 8.0115 × 10⁻⁹ m/s² between the results of the simplified and the simplified– modified algorithms. An improvement of about 1.8 m in the positions and 3 mm/s in the velocities is observed.

The MATLAB scripts developed in the context of this study were executed on an HP portable PC equipped with an Intel (R) Core (TM) i3 processor and installed memory (RAM) of 6.00 GB. The execution time (run time) of the simplified algorithm (six sub-equations) is 117.67 s and the simplified–modified algorithm (nine sub-equations) is 3778.33 s on average. The increase in the execution time is due to the addition of three equations and three initial conditions in the simplified–modified algorithm, and can be explained by the third figure (3d-diff A-L–S) in Figure 5. This sub-figure shows that the simplified–modified algorithm determines three new values of L–S accelerations every 30 min, differing from the values given in the broadcast ephemerides file, while the simplified algorithm uses directly the L–S accelerations as given in the broadcast ephemerides file.

Figure 5.

Difference between (Δ₁) and (Δ₂) of integrate orbits

In addition, the simplified–modified algorithm calculates every 30 min three linear functions representing the variation of the L–S acceleration along the three axes, whereas the simplified algorithm does not need these functions.

The use of the ode45 MATLAB solver has reduced the execution time compared to the execution time of another order, such as the 14th order, and this is one of the advantages of using this solver.

In the satellite-positioning concept, the unknown user’s position is determined from the known ephemeris of at least four satellites. The ephemeris of GLONASS satellites is not given at any desired time, but it is provided every 30 min, and for this, it is indispensable to compute these ephemerides at any desired time.

This paper analyzes the effect of the GLONASS satellites acceleration caused by the sun and moon perturbation in the computation of GLONASS orbits and exploits the advantages of the ode45 MATLAB solver based on the R–K orders 04 and 05 in the integration of the differential equation of GLONASS satellites’ motion.

To get the positions, velocities, and the L–S acceleration of GLONASS satellites, a simplified– modified algorithm was implemented in the ode45 MATLAB solver. The simplified algorithm given in the ICD-GLONASS (2016) document, which uses the L–S accelerations as a constant value during an interval of 15 min, has been modified to take these accelerations as a variable value during an interval of 30 min using the linear functions defined along the three axes. These functions give the derivatives of the L–S accelerations induced by the sun and moon perturbations on GLONASS satellites in orbits.

The implementation of the simplified–modified algorithm in the ode45 solver for GLONASS orbits computation eliminated the programming step of the K-coefficients given by equations 3 and 4 of the R–K method.

The integration errors have been diminished when the L–S accelerations are added in the simplified–modified algorithm as variable terms. The improvement is about 1.8 m in the position and 0.003 m/s in the velocity.

Another factor to take into consideration in the integration of GLONASS orbits by the R–K method is the computational time (run time). The ode45 solver minimizes this factor in a very important way since it is a built-in MATLAB function, and it is an advantage to use this solver in this study.