G i '• y GEODETSKI V E S T NIK | letu / Vol. 60 | št. / No. 3 V J 6n/3 3 RAZLIČNI PRISTOPI ZA DIFFERENT APPROACHES IZRAČUN TIRNIC GLONASS- IN GLONASS ORBIT SATELITOV IZ ODDANIH COMPUTATION FROM EFEMERID BROADCAST EPHEMERIS Kamil Maciuk UDK: 52-323.8 Klasifikacija prispevka po COBISS.SI: 1.01 Prispelo: 15. 3. 2016 Sprejeto: 12. 7. 2016 DOI: 10.15292/geodetski-vestnik.2016.03.455-466 SCIENTIFIC ARTICLE Received: 15. 3. 2016 Accepted: 12. 7. 2016 _ IZVLEČEK Enačbo gibanja umetnega satelita lahko rešujemo z različnimi metodami. Pri nekaterih diskretne položaje satelitov dobimo v enem, pri drugih v več korakih. Pri numerični integraciji na natančnost izračuna vpliva izbor velikosti koraka integracije med zaporednimi iteracijami. Nepravilen izbor integracijskega koraka med posameznimi iteracijami vodi do odstopanj, ki so lahko večja od dosegljive natančnosti izračuna s posamezno metodo. Zaradi tega moramo v izračunih integracijski korak toliko zmanjšati, da ne vpliva več na natančnost izračuna. V prispevku smo analizirali uporabnost različnih metod Runge-Kutta za numerično reševanje gibanja umetnega satelita, in sicer metodo Runge-Kutta 4. in 5. stopnje ter metodo Runge-Kutta-Fehlerberg 4. in 5. stopnje. Rezultate izračunov smo primerjali s klasično metodo Runge-Kutta in ugotovili, da z ustreznim izborom integracijskega koraka in metode nekoliko upočasnimo računski čas, vendar pridobimo bolj kakovostne rezultate izračunov. ABSTRACT _ Several types of methods can solve equations of satellite motion numerically. These methods are divided into single and multi-step methods. The accuracy of each method depends directly on adopted integration step size between successive iterations. To achieve result with required accuracy it is important to maintain appropriate size of integration step. Inappropriate step size could cause local errors between iterations greater than accuracy of the method. Therefore, integration step size needs to be reduced until it does not affect accuracy of the final solution. Group of Runge-Kutta (RK) methods for solving equations of satellite motion have been analysed in this article. Five different methods: Runge Kutta 4th order, Runge Kutta 5th order and Runge Kutta Fehlberg 4th and 5th order methods were discussed. Compared to the classical Runge-Kutta integration method other methods are slower, but give results that are slightly more accurate. KLJUČNE BESEDE KEY WORDS numerična integracija, GNSS, GLONASS, satelitska geodezija numerical integration, GNSS, GLONASS, satellite geodesy Kamil Maciuk | RAZLIČNI PRISTOPI ZA IZRAČUN TIRNIC GLONASS-SATELITOV IZ ODDANIH EFEMERID | DIFFERENT APPROACHES IN GLONASS ORBIT COMPUTATION FROM BROADCAST EPHE- MERIS |455-466 | I 455 | |60/3| GEODETSKI VESTNIK ^ 1 INTRODUCTION g Equations of satellite motion could be solved both analytically (Goral and Skorupa, 2012) and numeric ally (Gaglione et al., 2011). Runge-Kutta (RK) methods are one of the well-known numerical methods s for solving differential equations (Kosti et al., 2009; Ozawa, 1999; Sermutlu, 2004), while 4th order g Runge-Kutta method is recommended to solve equations of satellite motion by GLONASS Interface 1 Control Document (ICD-GLONASS, 2008). ^ Currently there are very few publications referring to comparison of numerical methods to solve GNSS equations of satellite motion. Numerical integration of low Earth orbiting satellites was performed by (Es-hagh, 2005). The author compared two variable step integration methods: Adams and RungeKut-taFehlberg (RKF) methods. Adams' method is recommended for long arc orbit integration or in low resolutions (large step size) orbit integration. In contrast, RKF method is better to be used for highresolution (small step size) solutions. (Sermutlu, 2004) presented the comparison of accuracy and speed tests of RungeKutta 4th and 5th order for solving Lorenz equation. He noticed that 4th order method gives more accurate results for shorter running times, but as step sizes decline, 5th order method gives more accurate results. (Khodabin and Rostami, 2015) obtained the same results. The authors analysed different orders of Runge-Kutta methods for applications in electric circuits. They confirmed superiority of higher order RK methods over other methods. (Montenbruck, 1992) compared multistep, interpolation and Runge-Kutta methods for the numerical integration of ordinary differential equations of orbital motion. The author showed that both single-step and multi-step methods are competitive. Equations of satellite motion were also solved by many different approaches, e.g. RungeKuttaFehlberg method (Atanassov, 2010), analytically (Kudryavtsev, 1995), by MATLAB ODE45 function (Bradley et al., 2014) or by new types of Runge Kutta methods (Gonzalez et al., 1999). Runge-Kutta 4th order method to solve equations of satellite motion was presented by (ICD-GLONASS, 2008), but without any data concerns accuracy. It is clear that the error in orbit integration strongly depends on a step size. GLONASS satellite integration results have no explicit differences between solutions from 1 to 300 s integration step size. The author suggested that 60 s GLONASS integration step width is sufficient in any case, because for small angular distances the satellite orbit could be considered as nearly linear. 2 KEPLERIAN MOTION Simplified satellite orbiting is called Keplerian motion (Zare, 1982). In Earth-artificial satellite, setting the mass of a satellite can be considered negligible and does not enter the motion equations system (Breiter and Elipe, 2006). This is due to its size and mass that are negligibly small relatively to the mass of the Earth. The satellites motion is governed by the Newton's second law hence, according to the formula: r=-4 r (i) where: fi= GM the product of Newton's gravitational constant and mass of the Earth, distance between the Earth and satellite centres. Kamil Maciuk | RAZLIČNI PRISTOPI ZA | 456 | MERIS | 455-466 | TIRNIC GLONASS-SATELITOV ) EFEMERID | DIFFERENT APPROACHES IN GLONASS ORBIT COMPUTATIC EPHE- r GEODETSKIVESTNIK | 60/3 | Equation 1 relates to a motion in an inertial system. Two vectors or 6 scalars are the solutions of this second order differential equation (Keplerian elements). They are the results of double integration of (1). In case of the Earth's artificial satellite, perturbing forces affecting its position should also be taken into account (Bobojc and Drozyner, 2011): - M r r = - —— + K (2) where: K perturbing forces. Gravitational forces due to the Earth as well as the strength of perturbing forces determine satellites motion. Table 1 shows the magnitude of perturbing forces and their effect on a GNSS satellite. Table 1: Perturbing accelerations acting on a GPS satellite (Dach et al., 2007). Source Acceleration [m/s2] Orbit error after 24 hours [m] Two-body term of the Earth's gravity field 0.59 ^ Oblateness of the Earth 5 • 10-5 10.000 Lunar gravitational attraction 5 • 10-6 3.000 Solar gravitational attraction 2 • 10-6 800 Other terms of the Earth's gravity field 3 • 10-7 200 Radiation pressure (direct) 9 • 10-8 200 Y-bias 5 • 10-10 2 Solid Earth tides 1 • 10-9 0.3 The main perturbing force affecting a satellite is the Earth's oblateness that characterizes polar flattening of the Earth. The effect of accelerations due to the luni-solar gravitational perturbations is an order of magnitude smaller than the second zonal harmonic. We can consider other forces as negligible. It may be assumed that perturbing forces acting on a GPS satellite affect will be different that on a GLONASS satellites due to two reasons. Firstly, GLONASS satellite orbit the Earth much lower, that is mean they are much sensitive to gravitational perturbations. Secondly, GLONASS satellites have larger area-to-mass ratios than GPS satellites, which implies that the impact of solar radiation pressure is larger for GLONASS. 3 RUNGE-KUTTA METHODS Numerical integration methods can be divided into single and multi-step methods. In case of multi-step methods, to calculate the predicted value of the function, we must know values of the function at some previous time points (e.g. tn , t 2). The best known multi-step methods used to solve equations of satellite motion are Cowell and Encke methods (Liu and Liao, 1994). Whereas single-step methods based on a single initial point of time, allow us to calculate predicted values of the function. The best-known singlestep methods for solving satellite equations of motion are Runge-Kutta 4th and higher order methods. The equation of satellite's motion is a second order differential equation. Therefore, it has to be converted to the system of first order differential equations to be solved by RK methods as following: Kamil Maciuk | RAZLICNIPRISTOPI ZA IZRACUN TIRNIC GLONASS-SATELITOV IZ ODDANIH EFEMERID | DIFFERENT APPROACHES IN GLONASS ORBIT COMPUTATION FROM BROADCAST EPHE- MERIS | 455-466 | | 457 |60/3| GEODETSKI VESTNIK y'(x) = f (x, y(x)) y(t0) = y (3) Runge-Kutta method allows calculation of the approximate value of the function y(xj for a = x„ < x. < ... < x = b, as in the formula: 0 1 n y , = y + h y b.k. Jn+1 Jn / • i i i=1 x , = x + h n+1 n (4) where: k,=f x +c h, y +h y a.k. n i ' J n / , i) ) and (5) )=! y = 1. ' (6) constants, step size. Runge-Kutta method's order, j=1 where: a, b -h -s - i=1,2,... ,s. Expanding (2) into first order differential equations still makes it impossible to solve them analytically in a fast and simple way. GLONASS Interface Control Document (ICD-GLONASS, 2008) recommends the use of Runge-Kutta 4th order method for this purpose, as it ensures adequate accuracy altogether with the simplicity of the solution. Equation (7) is an extension of (2) into a form of scalar functions. It takes into account perturbing forces due to the flattening of the Earth (second zonal harmonic) and influence of the Sun and Moon (Poutanen et al., 1996): dx — = x dt dy dt dz — = z dt dx Li 3 „ -= -—rx + — C. dt dy dt Li, 5z . .. ^ 0 20 5 x I 1--j" \ + xLS+a x + 2mJ 2 r I r (7) — = - — y + -C. r 3 T 20 5 r ua~ I 5z2 y|1-— r -yLs- _ 2a x dz, Li 3 „ — =--rZ + — C. dt 20 5 2r La2zi3 5z2 r r Kamil Maciuk | RAZLIČNI PRISTOPI ZA IZRAČUN TIRNIC GLONASS-SATELITOV IZ ODDANIH EFEMERID | DIFFERENT APPROACHES IN GLONASS ORBIT COMPUTATION FROM BROADCAST EPHE- 458 | MERIS | 455-466 | GEODETSKI VESTNIK | 60/3 | where: x, y, z x,y,z xLS ' JLS ' a a C„n satellite coordinates, satellite velocities, lunisolar accelerations, semi-major axis of the ellipsoid, the Earth rotation rate, second zonal harmonic coefficient of the geopotential, ■=4* y /u= GM. Second zonal harmonic is known from parameters of current PZ-90 (napaMeTpbi 3eMAH 1990 roga, Parameters of the Earth 1990) realization. In calculations, it is adopted as the known parameter. Lunisolar accelerations are varying in time, thus they are transmitted in GLONASS navigational (broadcast) message in 15 min intervals, and they are assumed constant within ±15 min from the initial position. LS 4 GLONASS NAVIGATION MESSAGE GLONASS navigation message contains information regarding satellites' position parameters for a single observation epoch. Those data are recorded in RINEX format *.yyG (Gurtner and Estey, 2007) with 30-minutes interval as vector components of satellite position, velocity and acceleration (Table 2). Table 2: GLONASS data record description (Gurtner and Estey, 2007). Observation record Description Format SV / EPOCH / SV CLK - Satellite system (R), satellite number A1, I2.2 (slot number in sat. constellation) - Epoch: Toc - Time of Clock (UTC) _- year (4 digits) 1X, I4 __ - month, day, hour, minute, second 5 (1X, I2, 2), - SV clock bias (sec) (-TauN) 3D19.12 - SV relative frequency bias (+GammaN) - Message frame time (tk+nd*86400) _ in seconds of the UTC week BROADCAST ORBIT - 1 - Satellite position X___(km) Satellit Satellit Satellit te velocity X dot__ te X acceleration_ te health (0=OK) 4X, 4D19.12 (km/sec) (km/sec2) BROADCAST ORBIT - 2 - Satellite position Y (km) 4X, 4D19.12 - Satellite velocity Y dot__ (km/sec) - Satellite Y acceleration__ (km/sec2) - Satellite frequency number_ _(-7... + 12) BROADCAST ORBIT - 3 - Satellite position Z _ (km) 4X, 4D19.12 - Satellite velocity Z dot__ (km/sec) - Satellite Z acceleration__ (km/sec2) - Age of oper. information__ _(days) Kamil Maciuk | RAZLIČNI PRISTOPI ZA MERIS | 455-466 TIRNIC GLONASS-SATELITOV ) EFEMERID | DIFFERENT APPROACHES IN GLONASS ORBIT COMPUTATIC EPHE- I 459 | |60/3| GEODETSKI VESTNIK Figure 1 contains a single record of GLONASS navigational message in RINEX format. It relates to the satellite PRN 1 from 9th June 2013 at 0:00 GLONASS time. PRN y m d h m s 1 13 6 9 0 0 0.0 Satellite position X (km) 0.144409179688E+05 Satellite position Y (km) 0.522635791016E+04 Satellite position Z (km) 0.203675307617E+05 SV clock bias (sec) -0.172111205757E-03 X velocity (km/sec) -0.264622497559E+01 Y velocity (km/sec) 0.877996444702E+00 Z velocity (km/sec) 0.165256118774E+01 SV relative frequency bias Message frame time 0.000000000000E+00 X acceleration (km/sec2) 0.000000000000E+00 Y acceleration (km/sec2) 0.000000000000E+00 Z acceleration (km/sec2) -0.279396772385E-08 0.846000000000E+05 Health 0.000000000000E+00 Frequency number 0.100000000000E+01 Age of oper. information 0.000000000000E+00 Figure 1: Example of GLONASS navigation message. Contrary to GPS, GLONASS message contains information about satellites' positions in ECEF coordinate system (Gaglione et al., 2011). Those data for a single satellite are stored in four 80-byte lines (Figure 1). GLONASS ephemeris message contains information about satellites' position in current PZ-90 realization (Boucher and Altamimi, 2001). PZ-90.02 realization was obligatory since 2007 (Montenbruck et al., 2015), currently PZ-90.11 is in use (IGSMAIL-6896). 5 GENERAL COMPARISON In this paper, a group of Runge-Kutta methods were analysed in resolving equations of satellite motion for GLONASS satellite. Parameters of GLONASS space segment are presented in Table 3. Table 3: GLONASS space segment parameters (Angrisano et al., 2013). Parameter Value Number of SV 24 Orbital planes 3 Orbital altitude (km) 19 100 Orbital inclination 64.8° Ground track period 8 sidereal days Layout Symmetric Broadcast ephemerides ECEF Datum PZ-90 This paper discusses four variants of Runge-Kutta method: best-known 4th order method (RK4), 5th order method (RK5) and Runge-Kutta-Fehlberg 4th (RKF4) and 5th (RKF5) order methods. Table 4 shows formulas of analysed RK methods. The determination error of satellite position depends on Runge-Kutta method order and adopted for calculations integration step. In principle, position determination is more accurate for smaller integration steps. Smaller integration step (h) carries a serious increase of intermediate positions thus, it increases computation time. Each step h, depending on the adopted formula requires calculation of four, five or six intermediate values of the function analysed in this paper (Table 4). Therefore, the best solution appears to be a method, which provides required accuracy of a satellite position solution combined with the highest execution speed. It is especially important in case of real-time solutions. | 460 | Kamil Maciuk | RAZLIČNI PRISTOPI ZA MERIS | 455-466 TIRNIC GLONASS-SATELITOV ) EFEMERID | DIFFERENT APPROACHES IN GLONASS ORBIT COMPUTATIC EPHE- GEODETSKI VESTNIK | 60/3 | Table 4: Parameter of analysed Runge-Kutta methods (Rentrop et al. 1989; Sermutlu, 2004). Runge-Kutta 4th order (RK4) Runge-Kutta 5th order (RK5) Runge-Kutta-Fehlberg (RKF45) k = hf (t„, J„ k = hf \ t + — h, y + — k 2 ^ n 2 n 2 1 k = hf \ K + 1 h, J, + 2 K k4 = hf (tn + h, Jn + k3 ) k + 24 + 24 + 4 Jn+1 = Jn + - 6 k, = hf (tn , Jn k2 = hf| K + 2h, Jn + 2k —6 k3 = hf I tn + 4 h, Jn = hf\ tn + 2 h, Jn + — k , 3, -3k2 + 6k, + 9k k5 = hf I tn + 4 h, Jn 2 ' 4 k6 = hf I tn + h, J, —6 4k + 6k -12k, Jn+1 = Jn 7 - 32k3 + 12k4 + 32k5 + 7k6 90 k, = hf (tn, Jn k = hf | tn + 4 h, Jn + 4 k 4 = hfl ^ 3, 3 , 9 . -—h, j +--k, +--k 8 /n 32 1 32 ...i 193^ 7200, = hf| t +--h, y +--k.--k 1 n 13 2197 1 2197 k5 = hfl tn + h, ', Jn = hf| t + —, y-- 1 n 2 27 Jn+1 = Jn 43^ ^ , 3680 --4 - 8k, H--k 216 1 2 513 3544, -k 2565 3 25 1408 2197 —k h--k h--k4 • 216 1 2565 3 4104 4 16 6656 28561 -4 H--k H--k 135 1 12825 3 56430 7296 +-k 2197 845 1 -4104' 1859 +-k 4104 1 5 -1i k 40 - — k 50 -— k 55 6 RESULTS This paper shows research of GLONASS' satellite position determination by RK methods according to the integration step size and its effect on the accuracy and speed of solution. There has been analysed position of #10 GLONASS satellite (SV 717, orbit 2, launched 25/12/2006, active from 03/04/2007) due date 01/01/2012 at three different moments of time: 1015, 1045 and 1115 UTC. The survey is based on broadcast orbit coordinates taken with maximum available accuracy of 12 decimal digits (Figure 1). Accuracy analysis was performed based on ORBGEN results, which is a part of Bernese GPS Software 5.0 (Dach et al., 2007). Comparison of numerical solutions of (2) was carried out based on the author's own scripts implemented in Matlab R2010b®. They were run on Lenovo L420 computer equipped with Windows 7 Professional, with the Intel Core i5-2410M 2.30 GHz, 4.00 GB RAM. Based on known initial function values of position, velocity and acceleration it is possible to determine satellite's position for any moment within the range ±15 min (900 s). This time span comes from the fact that the GLONASS ephemeris is updated every 30 minutes. If ephemeris data are used in the range exceeding ±15min difference between calculated and actual position expected is grow rapidly every ±15 min (Figure 2). 0 |15 Figure 2: Discrepancy of forward and backward 15 - minute integration. Kamil Maciuk | RAZLIČNI PRISTOPI ZA IZRAČUN TIRNIC GLONASS-SATELITOV IZ ODDANIH EFEMERID | DIFFERENT APPROACHES IN GLONASS ORBIT COMPUTATIC MERIS | 455-466 | EPHE- l 461 | |60/3| GEODETSKI VESTNIK Figure 3 shows errors of XYZ components calculated based on initial satellite position by RK4 with integration step h = 30 s. After 30 minutes, the error of each component does not exceed 1 meter, after 60 minutes error is smaller than 5 meters, and after 4 hours exceeds value of several meters. Therefore, in an application of numerical methods for solving equations of satellite motion information on satellite position in the shortest possible time intervals is very important. 75 - 50 25 1= o 0 UJ -25 -50 -75 .........o..........dX dY -a- dZ 0 30 60 90 120 150 180 210 240 Time [min] Figure 3: Increase of satellite position error (RK4, h = 30 s) Figure 4 shows more detailed data presented on Figure 3. "Known" coordinate and speed components are at t = 0 s. At t = 900 s follows update of satellite ephemeris data and then should be used next "known" position coordinate (t = 1800 s for this figure) and solved backward. Therefore, the increase of XYZ components error magnitude due to the updated ephemeris parameters is clearly visible. o UJ I 462 | 0 200 400 600 800 1000 1200 1400 1600 1800 Time [s] Figure 4: Increase of satellite position error (RK4, h = 1 s). Kamil Maciuk | RAZLIČNI PRISTOPI ZA IZRAČUN TIRNIC GLONASS-SATELITOV IZ ODDANIH EFEMERID | DIFFERENT APPROACHES IN GLONASS ORBIT COMPUTATIC MERIS | 455-466 | EPHE- GEODETSKI VESTNIK |60/3| Figure 5 shows the difference between RK4 method with the step h = 1 s solution and the reference solution. The figure presents three consecutive "backward" and "forward" solutions within 900 s interval. At 900 s, 2700 s and 4500 s moment, satellite coordinates, velocity and acceleration values are known. Solutions of three analysed, successive time points have similar errors. The offset of each component is a result of its update. That is why the determination of a single satellite position should be done within ±900 s around the known position. Maximum error in X component is around -0.1 m, Y around 0.9 m, and Z component up to 0.1 m error. Consequently, maximum 3D position error is 0.15 m. Thus, this type of calculation can be considered as sufficient for GLONASS broadcast orbit determination, due to its accuracy of about several meters. 0.2 o LU 0.1 0.0 -0.1 -0.2 \\ N!^......... / a \ jgs ^^^ / / -A- / y* jr 1 /v f ...................dX ---dY — dZ -3D x: \ 1000 2000 3000 Time [t] 4000 5000 \s Figure 5: Example of three consecutive integration steps. Table 5: Average duration of positions calculation and percentage changes in relation to RK4 [ms]. Step size h [s] 1 3 5 10 20 30 90 180 300 900 Number of steps 900 300 180 90 45 30 10 5 3 1 RK4 4.816 1.605 0.962 0.480 0.241 0.160 0.054 0.027 0.016 0.006 100% 100% 100% 100% 100% 100% 100% 100% 100% 100% RK5 7.374 2.881 1.758 0.992 0.468 0.353 0.147 0.084 0.060 0.038 153% 180% 183% 207% 194% 221% 272% 311% 375% 633% RKF4 7.712 2.539 1.557 0.776 0.401 0.275 0.107 0.066 0.049 0.031 160% 158% 162% 162% 166% 172% 198% 244% 306% 517% RKF5 8.047 2.719 1.608 0.774 0.403 0.275 0.113 0.079 0.054 0.033 167% 169% 167% 161% 167% 172% 209% 293% 338% 550% Table 5 presents a comparison of average speed of satellite position determination. These values are means of 100 000 consecutive solutions of Runge-Kutta methods. It depends on the adopted integration step size h. Increased integration step size decreases time of position determination. For each integration step the most efficient is Runge-Kutta 4th order method (RK4), due to the least complexity. The other three Kamil Maciuk | RAZLIČNI PRISTOPI ZA MERIS | 455-466 TIRNIC GLONASS-SATELITOV ) EFEMERID | DIFFERENT APPROACHES IN GLONASS ORBIT COMPUTATIC EPHE- I 463 | 0 |60/3| GEODETSKI VESTNIK ^ methods depending on the step length are between 2 to 6 times slower than RK4 method. Despite of P the most complex equations RKF5 method is the second fastest after the RK4 method among analysed. RKM is the slowest method for each step size. Speed of calculation in this method is comparable to other ^ only for 1 and 3 s integration step sizes. 1 Figure 6 presents calculated errors of "forward" satellite position. It reveals the difference between the =e author's and model solution based on integration step size. In case of small step length, less than 180 == s, results are comparable for all tested methods and the maximum error does not exceed 0.15 m. This accuracy is sufficient for navigation purposes. Integration step [s] Figure 6: Runge-Kutta method determination error [mm]. With the increase of integration step length a distinct advantage of higher order RungeKutta methods may be observed. It is clearly visible for integration steps h = 300 s and h = 900 s. RKF method projects satellite's trajectory with 0.60 cm accuracy for a single, 900 s step. If you need to determine denser number between consecutive positions/coordinates (e.g. coordinates are available every 60 s, but you want to have coordinates every 1 s) all you have to do is decrease step-size to needed. Therefore, simplicity is the main advantage of using this method against GPS, where navigation message data contain Keple-rian elements, which must have analytical solution. Moreover, GLONASS navigation message contains Cartesian coordinates and velocities in current PZ-90 realization every 30 min, so it is much affordable data than Keplerian elements in GPS navigation message available every 2 hours. 7 CONCLUSIONS The accuracy of GLONASS satellite's position calculated numerically depends mostly on integration step size. The influence of applied RK method type and order is smaller. Short integration step allows a relatively high precision, but it involves extension of solution time. Error of calculated position from initial parameter (epoch) increases together with "distance" from known coordinates. This study confirmed that higher order RK methods are more accurate. This fact is more evident especially in large-size Kamil Maciuk | RAZLICNIPRISTOPI ZA IZRACUN TIRNIC GLONASS-SATELITOV IZ ODDANIH EFEMERID | DIFFERENT APPROACHES IN GLONASS ORBIT COMPUTATION FROM BROADCAST EPHE- 4641 MERIS | 455-466 | GEODETSKI VESTNIK |60/3| integration steps of RK computations. The previous studies showed that the 5th order method or modified RKF methods are more accurate than the RK4 recommended by the GLONASS-ICD. On the other hand, due to the simplicity of equations RK4 order method is the fastest of the all analysed methods. However, an argument of economical saving time was more important in the 90s, when PCs' computing power was less efficient smaller than today. Currently due the highest accuracy of analysed methods, the most suitable for calculation of GLONASS satellite position is RungeKuttaFehlberg 5th order method. ACKNOWLEDGEMENTS The paper is the result of research carried out within statutory research no 11.11.150.006 and grant no. 15.11.150.006 in the Department of Geometrics, AGH University of Science and Technology, Cracow, Poland. Literature and references: Angrísano, A., Gaglíone, S., Gíoía, C. (2013). Performance assessment of GPS/GLONASS single point positioning in an urban environment. Acta Geodaetíca et Geophysíca, 48 (2), 149-161. DOI: http://dx.doi.org/10.1007/s40328-012-0010-4 Atanassov, A. M. (2010). Integration of the equation of the artificial Earth's satellites motion with selection of Runge-Kutta-Fehlberg schemes of optimum precision order. Aerospace Research in Bulgaria, 17, 24-34. Boucher, C., Altamímí, Z. (2001). ITRS, PZ-90 and WGS 84: current realizations and the related transformation parameters. Journal of Geodesy, 75 (11), 613-619. DOI: http://dx.doi.org/10.1007/s001900100208 Bradley, B. K., Jones, B. A., Beylkin, G., Sandberg, K., Axelrad, P (2014). Bandlimited implicit Runge - Kutta integration for Astrodynamics. Celestial Mechanics and Dynamical Astronomy, 119 (2), 143-168. DOI: http://dx.doi.org/10.1007/s10569-014-9551-x Breíter, S., Elipe, A. (2006). Critical inclination in the main problem of a massive satellite. Celestial Mechanics and Dynamical Astronomy, 95 (1), 287-297. DOI: http://dx.doi.org/10.1007/9781402053252_17 Dach, R., Hugentobler, U., Frídez, P, Meindl, M. (2007). Bernese GPS Software Version 5.0. Bern, Switzerland. Es-hagh, M. (2005). Step variable numerical orbit integration of a low earth orbiting satellite. Journal of the Earth & Space Physics, 31 (1), 1-12. Gaglíone, S., Angrísano, A., Puglíano, G., Robustellí, U., Santamaría, R., Vultaggio, M. (2011). A stochastic sigma model for GLONASS satellite pseudorange. Applied Geomatics, 3 (1), 49-57. DOI: http://dx.doi.org/10.1007/s12518-011-0046-0 Gonzalez, A. B., Martín, P, Lopez, D. J. (1999). Behaviour of a new type of Runge-Kutta methods when integrating satellite orbits. Celestial Mechanics and Dynamical Astronomyal Astronomy, 75 (1), 29-38. DOI: http://dx.doi.org/10.1023/A:1008387322426 Góral, W., Skorupa, B. (2012). Determination of intermediate orbit and position of GLONASS satellites based on the generalized problem of two fixed centers. Acta Geodynamica et Geomaterialia, 9 (3), 283-290. Gurtner, W., Estey, L. (2007). RINEX The Receiver Independent Exchange Format. ICD-GLONASS. (2008). GLONASS Interface Control Document,Version 5.1 (2). Moscow. Khodabin, M., Rostami, M. (2015). Mean square numerical solution of stochastic differential equations by fourth order Runge-Kutta method and its application in the electric circuits with noise. Advances in Difference Equations, 62 (1). DOI: http://dx.doi.org/10.1186/s13662-015-0398-6 Kosti, A. A., Anastassi, Z. A., Simos, T. E. (2009). An optimized explicit Runge-Kutta method with increased phase-lag order for the numerical solution of the Schrodinger equation and related problems. Journal of Mathematical Chemistry, 47 (1), 315-330. DOI: http://dx.doi.org/10.1007/s10910-009-9571-z Kudryavtsev, S. M. (1995). The fifth-order analytical solution of the equations of motion of a satellite in orbit around a non-spherical planet. Celestial Mechanics and Dynamical Astronomy, 61 (3), 207-215. DOI: http://dx.doi.org/10.1007/BF00051893 Liu, L., Liao, X. (1994). Numerical calculations in the orbital determination of an artificial satellite for a long arc. Celestial Mechanics and Dynamical Astronomy, 59 (3), 221-235. DOI: http://dx.doi.org/10.1007/BF00692873 Montenbruck, O. (1992). Numerical integration methods for orbital motion. Celestial Mechanics and Dynamical Astronomy, 53 (1), 59-69. DOI: http://dx.doi.org/10.1007/BF00049361 Montenbruck, O., Steigenberger, P., Hauschild, A. (2015). Broadcast versus precise ephemerides: a multi-GNSS perspective. GPS Solutions, 19 (2), 321-333. DOI: http://dx.doi.org/10.1007/s10291-014-0390-8 Bobojc, A., Drozyner, A. (2011). GOCE Satellite Orbit in the Aspect of Selected Gravitational Perturbations. Acta Geophysica, 59 (2), 428-452. DOI: http://dx.doi.org/10.2478/s11600-010-0052-3 Ozawa, K. (1999). A four-stage implicit Runge-Kutta-Nystrom method with variable coefficients for solving periodic initial value problems. Japan Journal of Industrial and Applied Mathematics, 16, 25-46. DOI: http://dx.doi.org/10.1007/BF03167523 Rentrop, P, Roche, M., Steinebach, G. (1989). The application of Rosenbrock-Wanner type methods with stepsize control in differential-algebraic equations. Numerische Mathematik, 55 (5), 545-563. DOI: http://dx.doi.org/10.1007/BF01398915 Kamil Maciuk | RAZLIČNI PRISTOPI ZA IZRAČUN TIRNIC GLONASS-SATELITOV IZ ODDANIH EFEMERID | DIFFERENT APPROACHES IN GLONASS ORBIT COMPUTATION FROM BROADCAST EPHE- MERIS | 455-466 | | 465 | 60/3 | GEODETSKI VESTNIK ^ Poutanen, M.,Vermeer, M., Jaakko, M. (1996).The permanent tide in GPS positioning. Journal of Geodesy, 70 (8), 499-504. | DOI: http://dx.doi.org/10.1007/BF00863622 | Sermutlu, E. (2004). Comparison of Runge-Kutta methods of order 4 and 5 on Lorenz equation. Journal of Arts and Sciences Say, 61-69. Zare, K. (1982). Numerical stabilization of keplerian motion. Celestial mechanics, 26(4), 407-412. DOI: http://dx.doi.org/10.1007/BF01230420 IGSMAIL-6896. (2014). https://igscb.jpl.nasa.gov/pipermail/igsmail/2014/008086.html. Maciuk K. (2016). Different approaches in GLONASS orbit computation from broadcast ephemeris. Geodetski vestnik, 60 (3): 455-466. DOI: 10.15292/geodetski- vestnik.2016.03.455-466 Kamil Maciuk, Ph.D. AGH University of Science and Technology Faculty of Mining Surveying and Environmental Engineering al. Mickiewicza 30 30-059 Krakow, Poland E-mail: maciuk@agh.edu.pl Kamil Maciuk | RAZLIČNI PRISTOPI ZA IZRAČUN TIRNIC GLONASS-SATELITOV IZ ODDANIH EFEMERID | DIFFERENT APPROACHES IN GLONASS ORBIT COMPUTATION FROM BROADCAST EPHE- 466 | MERIS | 455-466 |