# Determination of intermediate orbit and position of GLONASS satellites based on the generalized problem of two fixed centers.

1. INTRODUCTION

The problem of two fixed centers is a special case of the restricted three--body problem. The two fixed centers problem is well known in classical celestial mechanics: two fixed centers with masses [m.sub.1] and [m.sub.2], attract some massless particle, moving in their field according to Newton's laws. The integrability of this problem was the first time proved by Euler, by means of the separation of variables. A review of publications on the problem of two fixed centers, including its generalizations and applications are presented in (Lukyanov et al., 2005). A good example how to implement the model of the generalized problem of two fixed centers for geodetic applications is given in (Aksenov, 1969).

It is well known that the force function in the generalized problem of two fixed centers approximates the potential of the earth's gravitation within the square of the earth flattening (Aksenov, 1969; Lukyanov et al., 2005). It is therefore appropriate to adopt intermediate orbits for earth satellites obtained by solving the generalized problem of two fixed centers (Demin, 1970; Aksenov, 1977). The solution of the generalized problem of two fixed centers in a geocentric coordinate system may be written in the form (Aksenov, 1977, p. 49) :

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (1)

The fundamental plane of this system is the equator and the principal axis x points towards the vernal equinox. The rectangular coordinates x, y, z are associated with the spheroidal oblate coordinates [xi], [eta], and w (Aksenov, 1977, p. 49; Escobal, 1965, p. 146). The oblate spheroidal coordinate system finds substantial use in the theory of satellite motion based on the problem of two fixed centers. The coordinates [xi], [eta], and w are orthogonal: that is, the coordinate surfaces are an ellipsoid ([xi] = const.), a hyperboloid ([eta] = const.) and a plane (w = const.), which intersect one another at right angles (Demin, 1970; Aksenov, 1977). The constants c and [sigma] are specified by the conditions (Aksenov et al., 1963).

c = [a.sub.e] x [square root of [J.sub.2] - [([J.sub.3]/[2 x [J.sub.2]]).sup.2]], (2)

[sigma] = [[J.sub.3]/[2 x [J.sub.2]]] x [1/[square root of [J.sub.2] - [([J.sub.3]/[2 x [J.sub.2]]).sup.2]] = [[J.sub.3]/[2 x [J.sub.2]] x [[a.sub.e]/c] (3)

where [J.sub.2], [J.sub.3] are the second and the third gravitational zonal harmonics of the Earth, respectively, and [a.sub.e] is the equatorial radius of the reference ellipsoid. Assuming WGS-84, [a.sub.e] = 6378137.0 m, [J.sub.2] = 1.082626684 x [10.sup.-3], and [j.sub.3] = -2.53265649 x [10.sup.-6] (NIMA Technical Report, TR8350.2, 2000), we obtain: c = 209.728939 km, and [sigma] = -0.035571596.

2. CALCULATION OF THE INTERMEDIATE ORBIT ELLIPTIC ELEMENTS IN THE GENERALIZED PROBLEM OF TWO FIXED CENTERS

The intermediate satellite orbit in the generalized problem of two fixed centers is most simply described by the elements: a, e, i, [omega], [OMEGA], M which for c = 0 and [sigma] = 0 coincide with the corresponding Keplerian elements (Aksenov, 1977).

The position and velocity components of GLONASS satellites, taken from the navigation message, are re-computed from ECEF (Earth Centered, Earth Fixed) Greenwich coordinate system (PZ-90.02) to an absolute coordinate system.

We assume that at the initial moment t = [t.sub.0], the values of the geocentric rectangular coordinates in an absolute system [x.sub.0], [y.sub.0], [z.sub.0] as well as of their derivatives with respect to time [[??].sub.0], [[??].sub.0], [[??].sub.0], are known. The presented algorithms are in general based on the book of Aksenov, 1977.

First, from the following formulae three constants: [[alpha].sub.1] (the energy integral), [[alpha].sub.2] (the area integral), and [[alpha].sub.3] are determined (Aksenov, 1977).

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (4)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

From formulas (1) we find that

[[xi].sup.2.sub.0] = [[r.sup.2.sub.0] - [c.sup.2]]/2 [1 + [square root of 1 + [[4c.sup.2] [([z.sub.0] - c[sigma]).sup.2]/[[r.sup.2.sub.0] - [c.sup.2]).sup.2]]], (5)

[[eta].sub.0] = [[z.sub.0] - c[sigma]]/[[xi].sub.0], tg [w.sub.0] = [y.sub.0]/[x.sub.0] (6)

From (4) we obtain the following constants [??], [??], [??].

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (7)

The elements [??], [??] and [??] = sini are analogous to the major semi-axis, the eccentricity and the sine inclination of the elliptical orbit, respectively. When c = [sigma] = 0, they coincide. The parameters a, e, s can be determined by the method of successive approximations at the level of small quantities of order [[epsilon].sup.4].

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (8)

where the small dimensionless parameter [[epsilon].sub.i] is defined as follows:

[[epsilon].sub.i] = c/[a.sub.i](1 - [e.sup.2.sub.i]), i = 0,1,2,3 ...

For GPS and GLONASS satellites that parameter is less then 1/100.

Substituting [epsilon] = 0 we find approximate values for a, e, s, which we then adopt as first approximation. The process of approximation ends in case the following conditions: [absolute value of [a.sub.i] + 1 - [a.sub.i]] < [10.sup.-7], [absolute value of [e.sub.i+1] - [e.sub.i]] < 10-10, [absolute value of [s.sub.i+1] - [s.sub.i]] < [10.sup.-10], are fulfilled.

Practically the third approximation will be sufficiently accurate.

3. DETERMINING THE ANGULAR ELEMENTS [[omega].sub.0], [[OMEGA].sub.0] AND [M.sub.0] OF A SATELLITE ORBIT

First, according to (Aksenov, 1977, p. 108, 109) the auxiliary variables [[psi].sub.0] and [[theta].sub.0] are computed according to the following formulae:

cos[[PSI].sub.0] = a(1 - e[bar.e]) - [[xi].sub.0]/[[xi].sub.0] e - a ([bar.e] - e), (9)

sin[[psi].sub.0] = ae(1 - [[bar.e].sup.2])[J.sub.0][[??].sub.0]/[[sigma].sub.2][[[[xi].sub.0][bar.e] - a([bar.e] - e)].sup.2][square root of 1 - [k.sup.2.sub.2](1 - [cos.sup.2][[psi].sub.0])],

in which

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (10)

For [[THETA].sub.0] we have

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (11)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (12)

With [epsilon] = 0, [sigma] = 0, the angle [[psi].sub.0] coincides with the true anomaly and [[THETA].sub.0] with the argument of latitude of a Keplerian orbit. Next we calculate the orbital elements [[omega].sub.0], [[OMEGA].sub.0], and [M.sub.0] according to the following formulae:

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (13)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (14)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (15)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (16)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (17)

Next, from the following formula we obtain the eccentric anomaly [E.sub.0].

tg [[E.sub.0]/2] = [square root of [1 - [bar.e]]/[1 + [bar.e]]] tg [[psi].sub.0]/2], (18)

Finally we determine mean anomaly [M.sub.0].

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (19)

where the coefficients [lambda], [[lambda].sub.1], [[lambda].sub.2], [[lambda]'.sub.1], [[lambda]'.sub.2] and [e.sup.*] are defined as follows (Aksenov, 1977, p. 90):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (20)

With [epsilon] = 0, [sigma] = 0, the element [[omega].sub.0] coincides with the argument of perigee, [[OMEGA].sub.0] is the longitude of the ascending node and [M.sub.0] is the mean anomaly of the Keplerian orbit at epoch [t.sub.0].

4. ALGORITHM OF CALCULATION OF A SATELLITE POSITION AT INSTANT t

The position of the satellite in space is defined by the six known intermediate orbital elements a, e, s = sini, [[omega].sub.0], [[OMEGA].sub.0] and [M.sub.0].

M = [n.sub.0](t - [t.sub.0]) + [M.sub.0], (21)

[n.sub.0] = [([-2[alpha].sub.1]).sup.3/2]/GM (22)

The parameters [psi], [omega] and E are calculated by the technique of successive approximations:

tg [[[psi].sub.i+1]/2] = [square root of [1 + [bar.e]]/[1 - [bar.e]] tg [[E.sub.i]/2], (24)

[[omega].sub.i+1] = v[[psi].sub.i+1] + 1 + [[omega].sub.0]. (25)

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII]

A reasonable initial value for [E.sub.0] is M The process of approximation ends when the following condition is fulfilled: [absolute value of [E.sub.i+1] -[E.sub.i]] <[10.sup.-10] Usually 3-4 iterations are sufficient.

The Cartesian coordinates are computed (at instant t) using the following formulae (Aksenov, 1977, p.93):

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (27)

where

[MATHEMATICAL EXPRESSION NOT REPRODUCIBLE IN ASCII] (28)

5. NUMERICAL TESTS

Knowing the satellite position and velocity at epoch [t.sub.0] (Table 1) from broadcast ephemeris the intermediate orbital parameters of the GLONASS satellite GLN20 are computed. Detailed numerical examples of these computations are given in Table 2. Firstly, the coordinates X, Y, Z and velocity vector components (Table 1) at epoch [t.sub.0] are transformed from the ECEF coordinate frame (PZ-90.02) to an absolute coordinate system [x.sub.0], [y.sub.0], [z.sub.0], [[??].sub.0], [[??].sub.0], [[??].sub.0], (Table 2, 1/6) using the formulae given in GLONASS ICD 2008. Table 2 contains all the values of the computed parameters in agreement with the formulae and designations given in section 2 and 3. Computed orbital parameters are highlighted in bold. The algorithm of the computation of satellite positions, described in section 4, is illustrated in detail in Table 3 and Table 4.

Moreover, on the base of broadcast ephemeris the intermediate orbital elements of satellite GLN5 are computed for various epochs between 1:45 to 5:15 (March 18, 2011) at 30 minutes interval. Based on these orbital elements forward (for 15 minutes) and backward (for 15 minutes) positions of GLN5 satellite are computed. The difference between forward and backward positions are given in (Fig. 1). In the next numerical example, on the base of known orbital elements of GLN5 satellite obtained at epoch [t.sub.i] forward (for 30 minutes, [t.sub.i] + 30) and backward (for 30 minutes, [t.sub.i] - 30) satellite positions are calculated, and then the results were compared with the corresponding coordinates given in the broadcast ephemeris. Differences between these computations are shown in (Fig. 2). Worth mentioning is that from 21 coordinate differences (as shown in Fig. 2) only two differences exceed 2 m.

[FIGURE 1 OMITTED]

[FIGURE 2 OMITTED]

6. CONCLUSIONS

The major aim of this paper was to illustrate a procedure for computing orbital intermediate elements from the GLONASS broadcast ephemeris. The paper concentrates on the practical issues of implementing the model of the generalized problem of two fixed centers for computing satellite positions illustrated by the example of GLONASS broadcast ephemeris. The GLONASS broadcast ephemeris are transmitted as a half-hourly satellite state vector, expressed in PZ-90.02 geocentric Cartesian coordinates and are nominally valid for 30 minutes.

According to GLONASS ICD 2008 the computation of the satelite positron at time t from a satellite state vector with reference time [t.sub.e] requires the numerical integration where the longest integration period recommended should be 15 min forward and backward from the current state vector.

The proposed analytical algorithm provides good results when calculating positions of the GLONASS satellites in the time-range of 30 minutes forward and 30 minutes backward from the given current state vector. Results presented in Figure 2, are more accurate then the mean square error of broadcast positions of GLONASS-M satellites, which is about 10 m (GLONASS ICD, 2008, Table 4.2). The analytical solution, obtained from the generalized problem of two fixed centers, demands less time for calculation than the method of numerical integration.

ACKNOWLEDGEMENT

This work has been partially supported by the project No. 11.11.150. 006, AGH--University of Science and Technology in Krakow.

The authors are grateful to anonymous reviewers for their constructive comments and suggestions.

REFERENCES

Aksenov, P., Grebenikov, A. and Demin, V.G.: 1963, The generalized problem of motion about two fixed centers and its application to the theory of artificial Earth satellites, translated in Soviet Astron.--AJ. 7, No. 2, American Institute of Physics, 276-283.

Aksenov, E.P.: 1969, Mechanical interpretation of the force function in the generalized problem of two fixed centers, Soviet Astron.--AJ. 12, No 4, American Institute of Physics, 681-685.

Aksenov, E.P.: 1977, Theory of motion artificial Earth's satellites (in Russian), Nauka Press, Moscow, pp.360.

Demin,V.G.: 1970, Motion of an artificial satellite in an eccentric gravitation field, translated and published as NASA Technical translation, TT F-579, Wshington D.C. (Translation of Dvizheniye Iskusstvennogo Sputnika v Netsentral'nom Pole Tyagoteniya, Nauka Press, Moscow, 1968).

Emelyanov, N.V.: 1992, The acting analytical theory of artificial Earth satellites, Astronomical & Astrophysical Transactions, 1, 119 - 127, Gordon and Breach Science Publishers S. A.

Escobal, P.R.: 1965, Methods of orbit determination, J. Willey & Sons, Inc.

GLONASS Interface Control Document (ICD), Edition 5.1, Moscow, 2008.

Lukyanov, L.G., Emelyanov, N.V. and Shirmin, G.I.: 2005, Generalized problem of two fixed centers or the Darboux-Gredeaks problem, Cosmic Research, 43, No. 3, Springer Verlag, 186-191.

NIMA Technical Report, TR8350.2, 2000.

Wladyslaw GORAL (1) * and Bogdan SKORUPA (2)

(1) Bronislaw Markiewicz State Higher School of Technology and Economics in Jaroslaw

(2) AGH University of Science and Technology, Department of Geomatics, Al. Mickiewicza 30, 30-059 Krakow, Poland

* Corresponding author's e-mail: wgik@agh.edu.pl

(Received February 2012, accepted July 2012)
Table 1 Positions, velocities and luni-solar accelerations of the
GLONASS satellite GLN20, in the PZ-90.02 reference frame, from
broadcast ephemeris file, 17 October 2011. Epoch [t.sub.0] and
[t.sub.1].

words [units]      [t.sub.0] =           [t.sub.1] =
[11.sup.h]            [12.sup.h]
[45.sup.m]            [15.sup.m]
[00.0.sup.s] UTC      [00.0.sup.s] UTC

X [km]          12317.93408200        10391.4926758
Y [km]          -2245.13232422        3032.69384766
Z [km]          22212.8173828         23096.2607422
[??] [km/s]        -1.25356674194        -0.864193916321
[??] [km/s]         2.77420043945         3.04328060150
[??] [km/s]         0.980820655823       -0.00542259216309
[[??].sub.LS]      -0.931322574616        0.0
[km/              x [10.sup.-9]
[s.sup.2]]
[[??].sub.LS]       0.0                   0.00
[km/
[s.sup.2]]
[[??].sub.LS]      -0.186264514923       -0.279396772385
[km/              x [10.sup.-8]         x [10.sup.-8]
[s.sup.2]]

Table 2 Calculation intermediate orbital elements of the GLONASS
satellite GLN20 at epoch 17 October 2011,
[t.sub.0] = [11.sup.h] [45.sup.m] [00.0.sup.s] UTC.

item   words [units]

1      [x.sub.0] [km]         11881.413366

2      [y.sub.0] [km]         -3950.207035

3      [z.sub.0] [km]         22212.817383
4      [[??].sub.0]              -0.564121
[km/s]
5      [[??].sub.0]               3.788976
[km/s]
6      [[??].sub.0]               0.980821
[km/s]
7      [V.sub.0.sup.2]           15.6365838761
[[km.sup.2]/
[s.sup.2]]
8      [r.sub.0.sup.2]    650512863.357182
[[km.sup.2]]
9      [[xi].sub.0.       650502263.311214
sup.2]
[[km.sup.2]]
10     [[eta].sub.0]              0.8712144653
11     [r'.sub.0]               124.3079598242
[[km.sup.2]/s]
12     [J.sub.0]          650535649.531330
[[km.sup.2]]
13     [Q.sub.0]           -4137728.33716742
[[km.sup.4]/
[s.sup.2]]
14     [[alpha].sub.1]           -7.8132471441
[[km.sup.2]/
[s.sup.2]]
15     [[alpha].        10167603454.3809
sub.2.sup.2]
[[km.sup.4]/
[s.sup.2]]
16     [[alpha].sub.3]        42789.9984002332
[[km.sup.2]/s]
17     [??] [km]              25507.988833
18     [[??].sup.2]              -0.000010606
19     [[??].sup.2]               0.8199198099

20     [epsilon] (i=0)            0.0082220048

21     a (i=1) [km]           25507.67852
22     [e.sup.2] (i=1)            0.0000015658
23     [s.sup.2] (i=1)            0.8199097563

24     [epsilon] (i=1)            0.0082222049

25     a (i=4) [km]           25507.678491

26     [e.sup.2] (i=4)            0.0000015675

27     [s.sup.2] (i=4)            0.819909756
28     [epsilon] (i=4)            0.0082222049
29     [bar.e]                    0.0012519478
30     [[sigma].sub.2]       100834.261470836
[[km.sup.2]/s]
31     [[??].sub.0]               0.0049314043
[km/s]
32     [k.sub.2.sup.2]            0.00000000009
33     [[psi].sub.0]              1.4863757692
34     [[sigma].sub.1]       100835.150279399
[[km.sup.2]/s]
35     [[??].sub.0]               0.0000382876
36     [k.sub.1.sup.2]            0.0000554971
37     d                         -0.0002648328
38     [gamma]]                  -0.0001871300
39     [[THETA].sub.0]            1.2946375300
40     [[omega].sub.0]            6.0914509486
41     v                         -0.0000050595
42     [w.sub.0]                  5.9622124440
43     [mu]                      -0.0000430857
44     [[mu].sub.1]              -0.0000000718
45     [[mu].sub.2]              -0.00000000001
46     [[mu].sub.3]               0.0
47     [[mu]'.sub.0]             -0.0000000076
48     [[omega].sup.              6.0914434283
49     [[OMEGA].sub.0]            4.9805303799
51     [lambda]                   0.0000000012
52     [[lambda].                 0.0000000009
sub.1]
53     [[lambda].                 0.0
sub.2]
54     [[lambda]'.                0.0000000009
sub.1]
55     [[lambda]'.               -0.0000138571
sub.2]
56     e *                        0.0012519867
58     [n.sub.0]                  0.0001549725
[[s.sup.-1]]

Table 3 Luni-solar accelerations in an absolute coordinate system
and orbital parameters of the GLONASS satellite GLN20 at epoch
[t.sub.0] and [t.sub.1], 17 October 2011.

item    words [units]       [t.sub.0] =          [t.sub.1] =
11h 45m 00.0s UTC    12h 15m 00.0s UTC

1      [[??].sub.LS0]       -9.2212106887 x       0.0
[km/[s.sup.2]       [10.sup.-10]
2      [[??].sub.LS0]        1.3059277289 x       0.0
[km/[s.sup.2]       [10.sup.-10]
3      [[??].sub.LS0]       -1.8626451492 x      -2.7939677239 x
[km/[s.sup.2]       [10.sup.-9]          [10.sup.-9]
4      a (i=4) [km]      25507.678491         25507.66778
5      [e.sup.2] (i=4)       0.0000015675         0.0000015649
6      [s.sup.2] (i=4)       0.8199097560         0.8199097528
7      [[omega].sub.0]       6.0914509486         6.0911021308
8      [[OMEGA].sub.0]       4.9805303799         4.9805303167

Table 4 Calculation of position of the GLONASS satellite GLN20
at epoch t = 12h 00m 00.0s, 17 October 2011.

item   words [units]   t = [t.sub.0] + 15m   t = [t.sub.1] - 15m

2      [psi]               1.6245990807          1.6249490022
3      [omega]             6.0914427289          6.0910939094
4      E (i=1)             1.6246029369          1.6249528586
5      [psi]               1.6258529142          1.6262017771
6      [omega]             6.0914427226          6.0910939031
7      E (i=4)             1.6246028200          1.6249527413
8      [xi] [km]       25509.396008          25509.395025
9      [THETA]             1.4341122025          1.4341122460
11     [rho] [km]      25516.951759          25516.950778
12     [rho]'[km]      25516.09048           25516.089499
13     [alpha]             0.4243704089          0.4243704127
14     [beta]             -0.0002247771         -0.0002247772
18     x [km]          11259.895951          11259.895712
19     y [km]           -512.795156           -512.794967
20     z [km]          22876.805241          22876.804074