New Applications of Elliptic Functions and Integrals in GPS Inter-Satellite Communications with Account of General Relativity Theory

Abstract

During the last 15–20 years, the experimental methods for autonomous navigation and inter-satellite links have been developing rapidly in order to ensure navigation control and data processing without commands from Earth stations. Inter-satellite links are related to relative ranging between the satellites from one constellation or different constellations and measuring the distances between them with the precision of at least 1 $\mathsf{\mu }$m micrometer (=${10}^{-6}$ m), which should account for the bending of the light (radio or laser) signals due to the action of the Earth’s gravitational field. Thus, the theoretical calculation of the propagation time of a signal should be described in the framework of general relativity theory and the s.c. null cone equation. This review paper summarizes the latest achievements in calculating the propagation time of a signal, emitted by a GPS satellite, moving along a plane elliptical orbit or a space-oriented orbit, described by the full set of six Kepler parameters. It has been proved that for the case of plane elliptical orbit, the propagation time is expressed by a sum of elliptic integrals of the first, the second and the third kind, while for the second case (assuming that only the true anomaly angle is the dynamical parameter), the propagation time is expressed by a sum of elliptic integrals of the second and of the fourth order. For both cases, it has been proved that the propagation time represents a real-valued expression and not an imaginary one, as it should be. For the typical parameters of a GPS orbit, numerical calculations for the first case give acceptable values of the propagation time and, especially, the Shapiro delay term of the order of nanoseconds, thus confirming that this is a propagation time for the signal and not for the time of motion of the satellite. Theoretical arguments, related to general relativity and differential geometry have also been presented in favor of this conclusion. A new analytical method has been developed for transforming an elliptic integral in the Legendre form into an integral in the Weierstrass form. Two different representations have been found, one of them based on the method of four-dimensional uniformization, exposed in the monograph of Whittaker and Watson. The result of this approach is a new formulae for the Weierstrass invariants, depending in a complicated manner on the modulus parameter q of the elliptic integral in the Legendre form.

1. Introduction

In the past 10–15, years the problem of GPS (GPS means Global Positioning System) satellite-ground station communications has been replaced by the problem of autonomous navigation and inter-satellite communications (ISC) (links), which has been mentioned yet in 2005 in the monograph [1]. Autonomous navigation means that satellites should have the capability to transmit data between them via inter-satellite cross-link ranging [2] and, thus, to ensure navigation control and data processing without commands from Earth stations in the course of six months. It is important that next generation space missions will attain sub-millimeter precision of measuring distances beyond ${10}^{-6}$ m (1 micrometer = 1 $\mathsf{\mu }$m = ${10}^{-6}$ m) by means of ultrashort femtosecond pulse lasers. The theoretical description of such measurements is inevitably related to general relativity theory (GRT). This means that the relative ranging and relative velocity model should account also for the bending of the transmitting path of the signal, which is significant for such large distances between the satellites due to the action of the gravitational field. The bending of the light (laser) is related to the important and fundamental physical fact that due to the action of gravitational field, the signal travels a greater distance in comparison with the straight path distance in the absence of a gravitational field.

In the book [3], autonomous navigation is defined as referring to “processes in which the spacecraft without the support from the ground-based TT&C (telemetry and tracking control) system for a long time, relying on its onboard devices, obtains all kinds of measurement data, determining the navigation parameters like orbit, time and attitude”. In this monograph, DORIS and PRARE navigation systems are determined as non-autonomous, because “DORIS system can determine the spacecraft’s orbit with high accuracy, but both need to exchange information with the ground stations”. Thus, the main aim of autonomous navigation is to reduce the dependance of spacecrafts on the ground TT&C network and, in this way, to enhance the capability of the system’s anti-jam and autonomous survivability. It is evident also that the determination of the satellite’s exact location is achieved through the reception, processing and transmitting of ranging signals between the different satellites.

Autonomous GPS navigation systems are necessary, in view of the proposal of researchers from the University of Texas, Aerospace Corporation, National Bureau of Standards, International Business Machines Corporation (IBM) and Rockwell Automation Inc., to monitor nuclear explosions, based on the GPS inter-satellite communication link.

The development of inter-satellite laser communication systems in space technologies that enable super-high-speed data transfers at rates greater than 1 Gbps is widely applied also in cube-/nano-satellite platforms such as CubeLCT, AeroCube-7B/C, CLICK, LINCS-A/B, SOCRATES and LaserCube [4]. However, in order to establish laser communications, the high performance of the arc-second level pointing system is required. This is a difficulty for the nanosatellite platforms. It can be supposed why there is such a difficulty—since the data-transmission rate of 1 Gbps between the two nanosatellites is at an inter-satellite range of 1000 km, the effects of curving the trajectory of the laser signal may be considerable. In the publication [5], concerning the optical links of small satellites, it has been admitted that the s.c. “pointing error” arises not only due to tracking sensors and mechanical vibrations, but also due to the base motion of the satellite. This fact, together with the ranging of signals and the pointing error, illustrates the assumption that it is not only the motion of two satellites, but the effect of curving of the signal trajectory, that is important. The paper [5] also asserts that for a GPS constellation with 24 satellites, there are a total of 8–16 links—they can be forward and backward links in the same orbit, but as well as lateral links between adjacent orbits. Moreover, the distance of GPS inter-satellite cross-links can reach 49,465 km.

The establishment of inter-satellite links (ISL) and measurement communications is of primary importance for the relative ranging and relative velocity determination between the satellites from one constellation or from different constellations such as the GPS system, the Russian GLONASS, the European Galileo and the China Beidou second-generation system; all of them are considered to be inter-operable with each other.

The theory of inter-satellite communications (ISC) is developed in the series of papers by S. Turyshev, V. Toth, M. Sazhin [6,7] and S. Turyshev, N. Yu and V. Toth [8]. The theory in these three paper concerns the space missions GRAIL (Gravity Recovery and Interior Laboratory) [6], GRACE- FOLLOW-ON (GRACE-FO—Gravity Recovery and Climate Experiment—Follow On) mission [7] and the Atomic Clock Ensemble in Space (ACES) experiment [8,9,10] on the International Space Station (ISS).

2. A Brief Review of the Shapiro Delay Formulae and the Aims of This Paper

2.1. Standard Shapiro Delay Formulae, Derived from the Null ConeEquation

The basic theoretical approach to find the propagation time of signals in the near-Earth space with account also of the gravitational field and in the framework of general relativity theory is based on the s.c. Shapiro delay formulae, proposed in 1964 [11] and frequently used also in VLBI (very large baseline interferometry) (see the review article by Sovers, Fanselow and Jacobs [12]).

(1)$T_{AB}={\displaystyle \frac{R_{AB}}{c}}+{\displaystyle \frac{2G{M}_E}{c^3}}ln\left({\displaystyle \frac{r_A+r_B+R_{AB}}{r_A+r_B-R_{AB}}}\right),$

where the coordinates of the emitting and of the receiving satellite are correspondingly $\mid x_A\left(t_A\right)\mid =r_A$ and $\mid x_B\left(t_B\right)\mid =r_B$, and $R_{AB}=$$\mid x_A\left(t_A\right)-x_B\left(t_B\right)\mid$ is the Euclidean distance between the signal-emitting satellite and the signal-receiving satellite. Note that the space points of emission $x_A\left(t_A\right)$ and arrival $x_B\left(t_B\right)$ of the signal are taken at the time $t_A$ of emission and $t_B$ of arrival correspondingly. In (1), $G_{\oplus }{M}_{\oplus }$ is the geocentric gravitational constant, $M_{\oplus }$ is the Earth mass and $G_{\oplus }$ is the gravitational constant; $t=$ TCG denotes the s.c. geocentric coordinate time (TCG). Further (and also in the above formulae), we shall denote by $T_{AB}=T_B-T_A$ the signal propagation time and, thus, we have replaced the times $t_A$ and $t_B$ by $T_A$ and $T_B$, to underline that these are times of emission and reception. The second term in formulae (1) is the Shapiro time delay term, accounting for the signal delay due to the curved space-time. An understandable derivation of the Shapiro delay formulae can be found in the PhD dissertations [13,14].

In order to calculate the propagation time $T_{AB}$ of the signal, it should be remembered that according to general relativity (see the book [15]), the emitter of the first satellite emits a signal propagating on the gravitational null cone after setting up to zero the infinitesimal metric element $d{s}^2=0$:

(2)$d{s}^2=0=g_{00}{c}^{2}d{T}^2+2g_{oj}cdTd{x}^j+g_{ij}d{x}^{i}d{x}^j.$

For the null cone Equation (2), the solution of this quadratic algebraic equation with respect to the differential $dT$ can be given as follows [16]:

$dT=\pm {\displaystyle \frac{1}{c}}{\displaystyle \frac{1}{\sqrt[.]{-g_{00}}}}\sqrt[.]{\left(g_{ij}+{\displaystyle \frac{g_{0i}{g}_{0j}}{-g_{00}}}\right)d{x}^{i}d{x}^j}+$

(3)$+{\displaystyle \frac{1}{c}}\left({\displaystyle \frac{g_{0j}}{-g_{00}}}\right)d{x}^j,$

where the metric tensor components are determined for an Earth reference system with an origin at the centre of the Earth.

Further, the propagation time should be found after (1) choosing the metric tensor components $g_{ij}$, $g_{0i}$ and $g_{00}$ in the given expression for the zero infinitesimal metric $d{s}^2=0$ (2), and (2) choosing the coordinates $x^i$ ($i,j=1,2,3$) in the metric, which, additionally, may depend on one, two or more parameters. Note that we have assumed that $x^0=1$. In the present investigation, the coordinates shall be chosen to be related with the satellite’s motion along the elliptic orbit. Usually, geostationary satellites are situated on nearly circular orbits. Satellites from the GPS constellation have a very small eccentricity of the orbit of the order $e=0.01$; for the Russian constellation GLONASS the eccentricity is $e\approx 0.02$. Let us remember that if a and b are, respectively, the large and the small semi-axis of the ellipse, then by definition, $e={\displaystyle \frac{\sqrt{a^2-b^2}}{a}}={\displaystyle \frac{c}{a}}$ ($0\le e<1$), where c is the distance between the foci of the ellipse. Thus, the eccentricity measures to what extent the ellipse deviates from the circle.

2.2. The Purpose of This Paper and the Relation to theParametrization of the Coordinates in the Metric Tensor

The main goal of this paper will be the exact analytical calculation of the propagation time of a signal by using a solution of the null cone Equation (2) and the parametrization of the space coordinates of the satellite for two different cases.

2.2.1. Case A—Plane Elliptical Coordinates of the Satellite

This is the most simple case, when the two-dimensional $(x,y)$ coordinates of the plane elliptical orbit of the satellite are parametrized in terms of the large semi-major axis a of the ellipse, the eccentricity e and the eccentric anomaly angle E. This angle is a dynamical parameter, related to the motion of the satellite on the orbit. The analytical expression for the parametrization of the $(x,y)$ coordinates is

(4)$x=a(cosE-e),y=a\sqrt[.]{1-e^2}sinE.$

These formulaes are well-known in all basic courses on celestial mechanics (see, especially, [17,18]. The time coordinate is chosen to be the celestial time $t_{cel}$, which is related to the eccentric anomaly through the Kepler equation

(5)$E-esinE=n(t_{cel}-t_p)=M.$

In this equation, e is the eccentricity of the orbit, $n=\sqrt[.]{{\displaystyle \frac{GM}{a^3}}}$ is the mean motion, M is the mean anomaly and $t_p$ is the time of perigee passage for the satellite. The mean anomaly M is an angular variable, which increases uniformly with time and changes by ${360}^0$ during one revolution. The mean motion n has the following geometrical meaning: This is the motion of the satellite along an elliptical orbit, projected onto an uniform motion along a circle with a radius equal to the large axis of the ellipse. Usually, M is defined with respect to some reference time—this is the time $t_p$ of perigee passage, where the perigee is the point of minimal distance from the foci of the ellipse (the Earth is presumed to be at the foci). Since the eccentric anomaly angle E will play an important role in the further calculation, let us also remember how this notion is defined: If, from a point on the ellipse, a perpendicular is drawn towards the large axis of the ellipse, then this perpendicular intersects a circle with a centre O at a point P. Then, the eccentric anomaly represents the angle between the joining line $OP$ and the semi-major axis (the line of perigee passage).

The first main purpose of this paper will be to solve Equation (3) for $dT$ as a differential equation with respect to the eccentric anomaly angle E, for the given parametrization (4) of the coordinates $(x,y)$ and for the null cone metric of the near-Earth in the standard form

(6)$d{s}^2=-c^2\left(1+{\displaystyle \frac{2V}{c^2}}\right){\left(dT\right)}^2+\left(1-{\displaystyle \frac{2V_{.}}{c^2}}\right)\left({\left(dx\right)}^2+{\left(dy\right)}^2+{\left(dz\right)}^2\right)=0,$

where $V={\displaystyle \frac{G_{\oplus }{M}_{\oplus }}{r}}$ is the standard gravitational potential of the Earth. No account is taken of any harmonics due to the spherical form of the Earth since the GPS orbits are situated at a distance more than 20,000 km (considered from the centre of the Earth). In fact, the null cone metric is the same, derived from the standard metric for the near-Earth space in basic monographs of gravity theory such as the one by Fock [19].

2.2.2. First Aim of the Paper—Expressing the Propagation Time inTerms of a Sum of Elliptic Integrals of Different Kinds (First, Second andThird) for the Case of a Signal Emission by a Satellite on a PlaneElliptical Orbit

For the first time in [20,21] and in the review paper [22], the solution for the propagation time T is found to be expressed as a sum of zero-order elliptic integrals of the first, second and of the third kind. Later, it will become clear what the difference is between elliptic integrals of different orders and kinds. The imaginary unit appears in these elliptic integrals, so it is important to prove that the integrals are real-valued quantities, since the propagation time T is a real and not an imaginary quantity. The physical approximation that is used for finding the solution is $\beta ={\displaystyle \frac{2V}{c^2}}\ll 1$, which, physically, is related to the assumption about a weak gravity field. This assumption is consistent with the experimental data and the calculation of the numerical quantity $\beta =$${\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^{2}a}}=0.334\times {10}^{-9}\ll 1$ by some authors [23] and also with the relativistic effects of the order of $c^{-2}$. Further in the text, this accuracy will be denoted as $O\left(c^{-2}\right)$. This is a consequence of the experimental fact that a GPS constellation has to keep time to an accuracy of about four nanoseconds per day (4 ns/d = $4\times {10}^{-9}$ s/d) [23], which is compatible with the fractional time stability of $5\times {10}^{-14}$, maintained by atomic clocks. Naturally, the fractional stability has been constantly improving in the past and present years with the perspective of becoming smaller and smaller.

2.2.3. Accuracy of the Relativistic Effects and the Possibility ofFinding the Solution of the Null Cone Equation

The last fact is mentioned due to the following interesting consistency between the experimentally measured accuracy of relativistic effects and the found solution of the null cone Equation (6) for the propagation time T in terms of the eccentric anomaly angle E. It will be proved in Section 3.4 that it is possible for this solution to be expressed analytically only under the assumption $\beta ={\displaystyle \frac{2V}{c^2}}\ll 1$. The necessity to implement such an assumption, physically expressing the fact of the weak gravitational field, is related also to the fulfillment of the equivalence principle only for weak gravitational fields and slow motions—an important topic, discussed in the monograph by Fock [19].

2.2.4. Second Aim of the Paper—Theoretical Justification of theApproach from the Viewpoint of General Relativity Theory and DifferentialGeometry

The second main purpose of the paper is to justify the theoretical approach and to provide the reasoning for choosing as parametrization variables in the null cone Equation (6) the coordinates (4), parametrizing the satellites trajectory along the orbit in terms only of three variables, $a,e$ and E, among which only the eccentric anomaly angle E is the dynamical parameter. The need to justify the theoretical approach will refer not only to the simplified case of the Kepler coordinates (4), but also to the next (second) case of a space-oriented orbit, when the disposition of the orbit in space will be described by six Kepler parameters. They will be determined in the next section. From a conceptual point of view, this is a very non-trivial problem, since the presence of the Kepler coordinates (4) in the null cone Equation (6) may raise the question of whether it is correct to use coordinates, not related to the trajectory of the signal, but to the motion of the satellite. Following such wrong logic, then, it might turn out that the calculated time T is not a propagation time of a signal, but the time of motion for the satellite. Fortunately, this line of reasoning is not correct, because the essential fact is the null cone Equation (6) and not the parametrization of the space coordinates $x, y$ or $x, y, z$. This will be confirmed by (1) a theorem in the monograph by Fock [19], which, further, will be formulated in detail [22], and (2) a theorem from differential geometry [24] with a meaning that is very close to the meaning in the monograph by Fock (cited also in the papers [22,25]).

In fact, here, the problem can be stated in a different way, since the propagation time T will be found as a solution from the null cone equation (remember also that the Shapiro delay formulae (1) is also derived from the null cone equation), but it should be derived from the null geodesics. So, the question is as follows: Should the null geodesic equation (which, as known, is difficult to solve exactly) be used or one may use instead only the null cone equation? However, the null cone equation is compatible with the light-like geodesic equation in the sense that each solution of the null cone equation is a solution also of the light-like geodesic equation according to the theorem in the book by Fock [19], so it represents a geodesic. But then if it is a light-like geodesic, according to another theorem in [26], it is the shortest path (length) from one point of the geodesic to another point on the same geodesic. This is consistent with the Fermat principle about the propagation of light, asserting that light travels between two points for a minimal time [27]. Since the light velocity is constant, light will travel along the shortest path between the two points. In the book by Mishtenko and Fomenko [24], the theorem is formulated in the sense that each extremal of the null cone functional

(7)$E\left(\gamma \left(s\right)\right)=\int \mid {\stackrel{.}{\gamma }}^2\left(s\right)\mid ds=\int g_{ij}\left(x\right){\stackrel{.}{x}}^i\left(s\right){\stackrel{.}{x}}^j\left(s\right)ds$

is also an extremal of the length functional

(8)$L\left(\gamma \left(s\right)\right)=\int \sqrt[.]{g_{ij}\left(x\right){\stackrel{.}{x}}^i\left(s\right){\stackrel{.}{x}}^j\left(s\right)}ds,$

but the statement in the inverse direction is not always fulfilled. This means that the extremals of the null cone functional constitute a subset of the extremals of the length functional—see also the arguments in [22,25]. The length functional of the type (8) can be obtained also after the integration of the null cone Equation (6) for the case when the space coordinates $x,y,z$ are parametrized by several variables (see the next section)—in differential geometry, well studied are the s.c. first and second fundamental forms [24,26]. In subsequent papers, this will also be studied in detail, but this paper will focus only on finding the solutions for the propagation time in terms of elliptic integrals.

2.2.5. Third Aim of the Paper—Finding the Solution for thePropagation Time T for the Case of Space-Oriented Orbits in Terms ofElliptic Integrals of Higher Order (Second and Fourth)

Space-oriented) orbits are well known from celestial mechanics (further, we shall use the coordinate representation from the book [18]) and in this case, the position of the orbit in space is characterized by six Kepler parameters $(f,a,e,\Omega ,I,\omega )$, which are, in fact, the orbital coordinates. The analogues of the coordinate transformation (4) for this case of space-oriented orbits are (see the cite [28] and also the known monographs [29,30,31,32])

(9)$x={\displaystyle \frac{a(1-e^2)}{1+ecosf}}\left[cos\Omega cos(\omega +f)-sin\Omega sin(\omega +f)cosI\right],$

(10)$y={\displaystyle \frac{a(1-e^2)}{1+ecosf}}\left[sin\Omega cos(\omega +f)+cos\Omega sin(\omega +f)cosI\right],$

(11)$z={\displaystyle \frac{a(1-e^2)}{1+ecosf}}sin(\omega +f)sinI,$

where $r={\displaystyle \frac{a(1-e^2)}{1+ecosf}}$ is the radius vector in the orbital plane, the angle $\Omega$ of the longitude of the right ascension of the ascending node is the angle between the line of nodes and the direction to the vernal equinox, the argument of perigee (periapsis) $\omega$ is the angle within the orbital plane from the ascending node to the perigee in the direction of the satellite motion $(0\le \omega \le {360}^0)$ and the angle I is the inclination of the orbit with respect to the equatorial plane.

From celestial mechanics, it is also known that there is a relation between the eccentric anomaly angle E and the true anomaly angle f

(12)$tan{\displaystyle \frac{f}{2}}=\sqrt[.]{{\displaystyle \frac{1+e}{1-e}}}.tan{\displaystyle \frac{E}{2}}⟹f=2arctan\left[\sqrt[.]{{\displaystyle \frac{1+e}{1-e}}}.tan{\displaystyle \frac{E}{2}}\right].$

2.2.6. Elliptic Integrals of the Second and of the Fourth Order with anIntegration over the True Anomaly Angle f

Finding the solutions for the propagation time T of the signal in terms of elliptic integrals of the second and the fourth order (as shall be explained, the order of the elliptic integral is related to the order of the polynomial in the nominator) is the third main goal of this paper. Initially, it was performed in the paper [21] and some results have been summarized (with some new small additions and interpretations) in the review paper [22]. The mathematical properties of the space-oriented transformations (9)–(11) enable the definition of these elliptic integrals in terms of the rapidly changing dynamical parameter during the motion of the satellite—the true anomaly angle f. In fact, this should be considered as a first step in expressing the propagation time in terms of all the Kepler parameters $(f,a,e,\Omega ,I,\omega )$, characterizing the transformations (9)–(11). This is a much more complicated problem and for the moment, it is not clear whether the integrals over the other five variables $(a,e,\Omega ,I,\omega )$ will be elliptic ones or not. Naturally, during the propagation of the signal from one satellite to another, the true anomaly angle will be the most rapidly changing and the other five parameters will change very slowly.

Further, the square of the radius-vector $x^2+y^2+z^2$ will not contain any of the angles $\Omega ,I,\omega$. In fact, it can easily be calculated from the expressions (9)–(11) for $x,y,z$ that $x^2+y^2+z^2=r^2$, where

(13)$r={\displaystyle \frac{a(1-e^2)}{1+ecosf}}$

and also

(14)$dr={\displaystyle \frac{a(1-e^2)esinf}{{\left(1+ecosf\right)}^2}}.$

Thus, by its construction, r is a radius vector in the orbital plane, depending only on f as a dynamical variable. This made possible the calculation of the elliptic integral for the propagation time T in terms only of the variable f.

There is one more physical reason for the propagation time to be expressed only through the dynamical variable f. As shall be shown further, finding the propagation time T as a function of the eccentric anomaly E for the plane elliptical orbit is physically justifiable because the eccentric anomaly E is related to the celestial time $t_{cel}$ through the Kepler Equation (5) $E-esinE=n(t_{cel}-t_p)=M$.

Celestial time $t_{cel}$ can be introduced also for the motion of a satellite along a space-oriented orbit. If v denotes the celestial velocity and $d{t}_{cel}$-the infinitesimal celestial time, then $vd{t}_{cel}$ can be written as

$vd{t}_{cel.}=d{t}_{cel}\sqrt{{\left({\displaystyle \frac{dx}{d{t}_{cel}}}\right)}^2+{\left({\displaystyle \frac{dy}{d{t}_{cel}}}\right)}^2+{\left({\displaystyle \frac{dz}{d{t}_{cel}}}\right)}^2}$

(15)$=\sqrt{{\left({\displaystyle \frac{dx}{df}}\right)}^2+{\left({\displaystyle \frac{dy}{df}}\right)}^2+{\left({\displaystyle \frac{dz}{df}}\right)}^2}df=v_{f}df,$

where

(16)$v_f=\sqrt{1+e^2+2ecosf}$

is the velocity, associated with the true anomaly angle f.

2.2.7. Fourth Aim of the Paper: Finding New Analytical Methods forExpressing the Two Weierstrass Invariants in the Elliptic Integrals of ZeroOrder and in the Legendre Form

Essentially, in view of the application of elliptic integrals in GPS inter-satellite communications and also in several other areas of physics (cosmology, nonlinear equations, etc.), in the paper [33], general types of elliptic integrals were considered:

$J_n^{\left(4\right)}\left(y\right)=\int {\displaystyle \frac{y^{n}dy}{\sqrt{a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4}}}\mathrm{and}$

(17)$J_n^{\left(3\right)}\left(x\right)=\int {\displaystyle \frac{x^{n}dx}{\sqrt{a{x}^3+b{x}^2+cx+d}}}.$

Some more general types of transformations,

(18)$x={\displaystyle \frac{a}{y^2}}+b,$

were applied, so that a and b represent complicated functions of the modulus of the elliptic integrals. The transformation (18) has been applied only with respect to a special type of elliptic integrals—the s.c. elliptic integrals of zero order (i.e., $n=0$)—and in the Legendre form. This is the left integral in (19). The advantage of such an approach is that there is no restriction for the modulus of the elliptic integral to be small. Also, the theory of elliptic functions uses the standard parametrization of an elliptic curve of the type $y^2=4x^3-g_{2}x-g_3$ in terms of the Weierstrass function and its derivative, meaning that the function $x=\rho \left(z\right)$ and $y={\rho }^{\prime }\left(z\right)$ satisfy the above equation.

In this paper, it will be proved that the propagation time T will be expressed by means of elliptic integrals of the types (with some coefficients in front)

(19)$\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}\mathrm{and}\int {\displaystyle \frac{y^{2}dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}.$

Since it is known that integrals in the Weierstrass form $\int {\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}$ can be solved analytically, it is evident that the transformation (18) can be used to find analytical expressions for the first integral in (19), which is called an elliptic integral of the first kind and of zero order in the Legendre form.

The second integral in (19) can also be solved analytically due to the recurrent chain of elliptic integrals of different orders, which usually is written as [34]

$y^{n}Z=a_0(n+2)J_{n+3}+2a_1(2n+3)J_{n+2}$

(20)$+6a_2(n+1)+2a_3(2n+1)J_n+n{a}_4{J}_{n-1},$

where Z denotes the square of the fourth degree polynomial in the integral $J_n^{\left(4\right)}\left(y\right)$ in (17).

(21)$Z:\equiv \sqrt{a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4}.$

Further, it will be shown that the invariant $g_2$ and the Weierstrass function $\rho (z,g_2,g_3)$ can be found from a complicated system of two nonlinear equations. Then, if the function $\rho (z,g_2,g_3)$ is known, one can apply the formulas (217)–(220), given in the monograph by Whittaker and Watson [35] with the purpose of finding a parametrization of the polynomials (21) and (201) $Z^2=f\left(y\right)=a_0{y}^4+a_1{y}^3+a_2{y}^2+a_{3}y+a_4=Y$ in terms of expressions, depending on the Weierstrass function and its derivative.

3. Propagation Time for the Case of a Signal-Emitting Satellite,Moving Along a Plane Elliptical Orbit

3.1. New Calculation of the Propagation Time in Terms of theEccentric Anomaly Angle E Without Any Approximations

As mentioned in the introduction of the paper, the null cone metric is taken in the standard form (6) (see also [36,37]):

$d{s}^2=-c^2\left(1+{\displaystyle \frac{2V}{c^2}}\right){\left(dT\right)}^2+\left(1-{\displaystyle \frac{2V_{.}}{c^2}}\right)\left({\left(dx\right)}^2+{\left(dy\right)}^2+{\left(dz\right)}^2\right)=0,$

where $V={\displaystyle \frac{G_{\oplus }{M}_{\oplus }}{r}}$ is the standard gravitational potential of the Earth. No account is taken of any harmonics due to the spherical form of the Earth since the GPS orbits are situated at a distance more than 20,000 km (considered from the centre of the Earth). The Kepler parametrization of the space coordinates $x=x(a,e,E)$ and $y=y(a,e,E)$ for the plane elliptical satellite orbit is (4) $x=a(cosE-e)$, $y=a\sqrt[.]{1-e^2}sinE$, where a is the semi-major axis of the orbit, e is the eccentricity and E is the eccentric anomaly angle, related to the motion of the satellite along the plane elliptical orbit.

If the right-hand side of Equation (6) is divided and multiplied by ${\left(dt\right)}^2$, where $t=t_{cel}$ is the celestial time from the Kepler equation $E-esinE=n(t_{cel}-t_p)$ ($t_p$ is the time of perigee passage), one can obtain

(22)$(c^2+2V){\left(dT\right)}^2=(1-{\displaystyle \frac{2V}{c^2}})v^2,$

where $v^2$ is the square of the satellite velocity along the orbit:

(23)$v^2=v_x^2+v_y^2={\left({\displaystyle \frac{dx}{dt}}\right)}^2+{\left({\displaystyle \frac{dy}{dt}}\right)}^2.$

Taking into account the Kepler parametrization (4), the above expression for the velocity for the case of plane motion can be rewritten as

(24)$v=\sqrt{v_x^2+v_y^2}={\displaystyle \frac{na}{\left(1-ecosE\right)}}\sqrt[.]{1-e^2{cos}^{2}E},$

where $n=\sqrt[.]{{\displaystyle \frac{G_{\oplus }{M}_{\oplus }}{a^3}}}$ defines the mean motion. Expressing $dT$ from (22), finding the celestial time $t_{cel}$ in terms of E from the Kepler equation and performing the integration over the eccentric anomaly angle E, one can obtain the expression for the propagation time (for the case of no restrictions imposed):

$T=\int {\displaystyle \frac{v}{c}}\sqrt[.]{{\displaystyle \frac{(c^2-2V)}{(c^2+2V)}}}dE+C$

(25)$={\displaystyle \frac{a}{c}}\int \sqrt[.]{{\displaystyle \frac{(1-e^2{cos}^{2}E)\left[a\left(1-ecosE\right)-\beta \right]}{\left[a\left(1-ecosE\right)+\beta \right]}}}dE+C.$

The constant $\beta$ is determined as

(26)$\beta ={\displaystyle \frac{2G{M}_{\oplus }}{a{c}^2}}.$

The other constant C can be determined as a result of the integration of the null cone equation from some initial condition—for example, the propagation time T should be equal to zero, when $E=0$ (if $E_{init}=0$). Further, making the substitution

(27)$1-ecosE=y\Rightarrow dE={\displaystyle \frac{dy}{esinE}}$

in the integral (25) and denoting $\overline{\beta }={\displaystyle \frac{\beta }{a}}$, one can represent (25) as

(28)$T={\displaystyle \frac{a}{c}}\int \sqrt[.]{{\displaystyle \frac{y(2-y)(y-\overline{\beta })}{(y+\overline{\beta })\left[e^2-{\left(1-y\right)}^2\right]}}}dy.$

This integral is of the form

(29)$T={\displaystyle \frac{a}{c}}\int \sqrt[.]{{\displaystyle \frac{a_1{y}^3+a_2{y}^2+a_{3}y}{b_1{y}^3+b_2{y}^2+b_{3}y+b_4}}}dy,$

where the numerical constants $a_1$, $a_2$, $a_3$, $b_1$, $b_2$, $b_3$, $b_4$ can easily be calculated. Their analytical expressions are not given here, because the above expression will not be used further.

3.2. Abelian Integrals Versus Elliptic Integrals

3.2.1. Definitions About Abelian Integrals

Abelian integrals are defined in the monograph [38] as

(30)$\int R(x,\sqrt[m]{{\displaystyle \frac{\alpha x+\beta }{\gamma x+\delta }}})dx,\int R(x,\sqrt{a{x}^2+bx+c})dx.$

The reason, according to the definition in [38], is that the integrals (30) are considered to be related to algebraic curves of the type

(31)$\left(\gamma x+\delta )y^m-(\alpha x+\beta \right)=0,$

(32)$y^2-(a{x}^2+bx+c)=0,$

obtained after setting up

(33)$y=\sqrt[m]{{\displaystyle \frac{\alpha x+\beta }{\gamma x+\delta }}}\mathrm{or}y=\sqrt{a{x}^2+bx+c}.$

It has been pointed out in [38] that abelian integrals of the type

(34)$\int {\displaystyle \frac{P\left(x\right)dx}{\sqrt[.]{a{x}^2+bx+c}}},\int {\displaystyle \frac{Adx}{{(x-\alpha )}^k\sqrt[.]{a{x}^2+bx+c}}},$

(35)$\int {\displaystyle \frac{\left(Mx+N\right)dx}{{(x^2+px+q)}^m\sqrt[.]{a{x}^2+bx+c}}}$

($P\left(x\right)$-polynomial, where $A,M,N$ are numerical constants) can be treated analytically, the second integral, for example, by means of the known Euler substitutions. Evidently, this is possible since the degree of the polynomial under the square root in the denominator of these integrals is two.

3.2.2. Is the Integral (29) an Abelian One?

If the square root, representing an irrational function of the variable y, is denoted as

(36)$x=\sqrt[.]{{\displaystyle \frac{a_1{y}^3+a_2{y}^2+a_{3}y}{b_1{y}^3+b_2{y}^2+b_{3}y+b_4}}}=\sqrt[.]{P\left(y\right)},$

then the integral

(37)$\int xdy=\int \sqrt[.]{P\left(y\right)}dy$

is a partial case of a class of integrals, often denoted in mathematical literature as $\int R(x,y)dy$). In case $P\left(y\right)$ is an algebraic polynomial, integrals $\int \sqrt[.]{P\left(y\right)}dy$ are known in mathematics as abelian integrals (see the monograph by Prasolov and Solovyev [34]), related to the curve

(38)$F(x,y):=x^2-P\left(y\right)=0.$

Since in the case of the expression (36), $P\left(y\right)$ is a complicated irrational function of an rational expression of y, the integral (29) is not an abelian one. Further, the expression Equation (25) will not be used, because a justifiable and reasonable physical approximation shall be applied, which will enable us to express the solution in terms of elliptic integrals. Nevertheless, in view of the constantly improving accuracy for measuring the propagation time, this might represent an interesting problem for future mathematical research. In the next sections, when defining the elliptic integrals for the more general case of arbitrary polynomials of the third and fourth degree in the denominator, it will be explained in detail why integrals of the form (29) cannot belong to the class of abelian integrals.

3.3. The Useful Approximation for a Small Gravitational Potential,Compared to the Square of the Velocity of Light

Further, we shall be interested in the case

(39)$\beta ={\displaystyle \frac{2V}{c^2}}={\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^{2}a}}\ll 1,$

which can be assumed to be fulfilled and is consistent with the current experimental data.It shall be proved in the next section that this physical assumption will greatly facilitate the calculations so that the appearance of any non-abelian integrals will be avoided.

For the parameters of the $GPS$ orbit—$a=26,561$ [km], in the review paper [23], the constant $\beta$ can exactly be calculated to be $0.334\times {10}^{-9}$, with the velocity of light taken to be $c=299,792,458$ $\left[{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}\right].$ The geocentric gravitational constant $G_{\oplus }{M}_{\oplus }$ (obtained from the analysis of laser distance measurements of artificial Earth satellites) can be taken to be equal to $G_{\oplus }{M}_{\oplus }=\left(3,986,004.405\pm 1\right)\times {10}^8\left[{\displaystyle \frac{{\mathrm{m}}^3}{{\mathrm{s}}^2}}\right]$, but the value of $G_{\oplus }{M}_{\oplus }$ can vary also in another range from $G_{\oplus }{M}_{\oplus }=3,986,056.75236\times {10}^8\left[{\displaystyle \frac{{\mathrm{m}}^3}{{\mathrm{s}}^2}}\right]$ to the value $G_{\oplus }{M}_{\oplus }=3,987,999.07898\times {10}^8\left[{\displaystyle \frac{{\mathrm{m}}^3}{{\mathrm{s}}^2}}\right]$ due to the uncertainties in measuring the Newton gravitational constant $G_{\oplus }$. The mass of the Earth can be taken approximately to be $M_E\approx 5.97\times {10}^{24}$ [kg]. One of the latest values for $G_{\oplus }$ from deep space experiments was reported in the paper [39] to be $(6.674+0.0003)$.${10}^{-11}$$\left[{\displaystyle \frac{{\mathrm{m}}^3}{\mathrm{kg}.{\mathrm{s}}^2}}\right]$. The latest value for G from Encyclopedia Britanica (updated 23 September 2024) is $6.67\pm 0.00015\times {10}^{-11}$$\left[{\displaystyle \frac{{\mathrm{m}}^3}{\mathrm{kg}.{\mathrm{s}}^2}}\right]$.

3.4. Derivation of the Expression for the Propagation Time Under theApproximation $\beta ={\displaystyle \frac{2V}{c^2}}\ll 1$ and PhysicalJustification of the Obtained Result

After decomposing the under-integral expression with the square root in Equation (29) under the assumption (39) and leaving only the first-order term in ${\displaystyle \frac{2V}{c^2}}$, one can obtain

(40)$T=\int {\displaystyle \frac{v}{c}}\sqrt[.]{{\displaystyle \frac{(1-\frac{2V}{c^2})}{(1+\frac{2V}{c^2})}}}dt\approx \int {\displaystyle \frac{v}{c}}(1-{\displaystyle \frac{2V}{c^2}})dt=I_1+I_2$

(41)$={\displaystyle \frac{a}{c}}\int \sqrt[.]{1-e^2{cos}^{2}E}dE-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}\int \sqrt[.]{{\displaystyle \frac{1+ecosE}{1-ecosE}}}dE.$

This expression is consistent from a physical point of view due to the following reasons:

The coefficient ${\displaystyle \frac{a}{c}}$ as a ratio of the large semi-major axis of the orbit and the velocity of light $c=299,792,458$ $\left[{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}\right]$ will have a dimension $[\mathrm{m}/{\displaystyle \frac{\mathrm{m}}{\mathrm{s}}}]=\left[\mathrm{s}\right]$, as it should be. The second coefficient ${\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}$ has a corresponding dimension $[{\displaystyle \frac{{\mathrm{m}}^3}{{\mathrm{s}}^2}}:{\displaystyle \frac{{\mathrm{m}}^3}{{\mathrm{s}}^3}}]=\left[\mathrm{s}\right]$, which clearly proves that formulae Equation (41) has the proper dimensions.

The formulae, in fact, gives the propagation time T of the signal, emitted by the satellite at some initial position (given by the eccentric anomaly angle $E_{init}$), and the final point (given by $E_{fin}$) on the trajectory of the signal. However, there is no real transmission of signals, because the final point $E_{fin}$ should be given. Let us clarify this moment—if there is a transmission of signals (emission+reception) the point of reception should be determined, depending on the equations for the two null cones in the framework of the approach of the two intersecting null cones, which was developed in the papers [20,40] for the simplest case of plane elliptic orbit. The method shall be applied for the case of space-oriented orbits in another paper. The integrals (40) and (41) represent curvilinear integrals (depending on the determined-in-advance initial and final points of integration), derived after the application of the null cone equation (in a general form, written as $F=g_{\alpha \beta }d{x}^{\alpha }d{x}^{\beta }=0$). This propagation time, in fact, can be found from the expression (3) after expressing the differentials $dT$ and $d{x}_i$ by means of the differential $dE$ of the eccentric anomaly angle E. Thus, from (3), after integration, it is obtained for the propagation time T:

$T=\pm {\displaystyle \frac{1}{c}}\int {\displaystyle \frac{1}{\sqrt[.]{-g_{00}}}}\sqrt[.]{\left(g_{ij}+{\displaystyle \frac{g_{0i}{g}_{0j}}{-g_{00}}}\right){\displaystyle \frac{d{x}^i}{dE}}{\displaystyle \frac{d{x}^j}{dE}}}dE$

(42)$+{\displaystyle \frac{1}{c}}\int \left({\displaystyle \frac{g_{0j}}{-g_{00}}}\right){\displaystyle \frac{d{x}^j}{dE}}dE.$

In the same way, the space variables $x^i$ can be parametrized by some other variable. For example, in the next section, by using the celestial transformations (9)–(11) for the case of space-oriented orbits, the true anomaly angle f will be used and the formulae for T will be the same as in (42) but with the eccentric anomaly E replaced by the true anomaly f, i.e., $E\to f$. This case will be investigated in the next sections. In order to obtain the propagation time T in the form of the expression (41), the additional assumption $\beta ={\displaystyle \frac{2V}{c^2}}\ll 1$ has to be taken into account, when transforming expression (42).

The expression for the propagation time T (41) establishes the relation between the eccentric anomaly angle E and the propagation time $T-$$E⟹T$. Since T will be expressed through elliptic integrals and their analytical solutions will be complicated, it might turn out to be impossible to establish the correspondence in the other direction $T\Rightarrow E$, as was in the case of the “dual” correspondence $t_{cel}\iff E$. The meaning of the correspondence $E⟹T$ is the following: If some initial and final eccentric anomaly angles $E_{init}$ and $E_{fin}$ for the position of the satellite are given, then the propagation time of the signal, corresponding to these two positions of the satellite, can be found from the definite integral

$\Delta T=T_{fin}-T_{init}={\displaystyle \frac{a}{c}}{\int }_{E_{init}}^{E_{fin}}\sqrt[.]{1-e^2{cos}^{2}E}dE$

(43)$-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}{\int }_{E_{init}}^{E_{fin}}\sqrt[.]{{\displaystyle \frac{1+ecosE}{1-ecosE}}}dE.$

This expression establishes the correspondence ($E_{init},E_{fin}$) $⟹T$, accounting for the curved signal trajectory and also for the fact that there is no real transmission of the signal, in the sense that the signal is not intercepted by another (second) satellite.

It should be stressed that T was found as a solution of the null cone Equation (6) and because of that, it is the time for propagation of the signal and definitely does not represent the time, related to the real motion of the satellite. This statement is also confirmed by the fact that the null cone equation $F=0$ is proved to be compatible with the geodesic equation [19] for light-like geodesics:

(44)${\displaystyle \frac{d^2{x}_{\nu }}{d{\tau }^2}}+{\Gamma }_{\alpha \beta }^{\nu }{\displaystyle \frac{d{x}^{\alpha }}{d\tau }}{\displaystyle \frac{d{x}^{\beta }}{d\tau }}=0,$

where $\tau$ is the proper time along the light ray and $\nu ,\alpha ,\beta =0,1,2,3$. In other words, the zero-length (light-like) geodesics (representing the trajectory of the signal) are determined also by the null cone equation $g_{\alpha \beta }d{x}^{\alpha }d{x}^{\beta }=0$, because the null cone equation is a first integral of the geodesic Equation (44). This point will be clarified in more detail in a separate chapter, since it is confirmed also by a theorem in differential geometry [24,26]. This was explained also in the introduction of this paper and also in the recent paper [25]. Furthermore, the numerical calculations in the next sections will confirm that the propagation time of the signal is much smaller than the time of motion of the satellite from the initial value $E_{init}$ of the eccentric anomaly angle to the final value $E_{fin}$.

Now, let $x,y,z$ be the space coordinates, reached by the signal after it has been emitted. From the null cone Equation (6), after using the approximation $\beta ={\displaystyle \frac{2V}{c^2}}\ll 1$, it can be obtained that

(45)$\left({\left(dx\right)}^2+{\left(dy\right)}^2+{\left(dz\right)}^2\right)\approx c^2\left[1+{\displaystyle \frac{4V}{c^2}}\right]{\left(dT\right)}^2\ll 3c^2{\left(dT\right)}^2.$

Since $cdT$ is the infinitesimal distance, travelled by the light signal, the above inequality means that the distance travelled by light is much greater than the Euclidean distance. In other words, the light trajectory signal is a curved one, in accord with the meaning of the Shapiro delay formulae (1) regarding the delay of the signal under the action of the gravitational field.

The trigonometric function ${cos}^{2}E$ has the same values for $E=\alpha$ and $E=180-\alpha$, so the correspondence eccentric anomaly angle E and the propagation time T are not reliable for angles in the second quadrant. The same refers to values of E in the first and the fourth quadrant since $cosE=cos(360-E)$.

The propagation time T in Equation (41) should be a real-valued expression, since time cannot be complex. This fact is not evident from the beginning of the calculations, but it will be proved in the next sections. Remarkably, this important property will turn out to be valid also for the next case of a signal, emitted and received by satellites on a space-oriented orbit.

The approximation $\beta ={\displaystyle \frac{2V}{c^2}}\ll 1$ with account of the vector r in the orbital plane is

(46)${\displaystyle \frac{2V}{c^2}}={\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^{2}a(1-ecosE)}}\ll 1,$

which is fulfilled, because $\beta =$${\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^{2}a}}\ll 1$. For the parameters of a $GPS$ orbit, the constant $\beta$ can exactly be calculated to be $0.33\times {10}^{-9}$, which justifies the above strong inequality. In the strict mathematical sense and for the given case, the above inequality (46) takes place when $cosE\ll {\displaystyle \frac{1}{e}}-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^{2}ae}}$. This is always fulfilled for the case of the small eccentricity of the GPS orbits ($e=0.01323881349526$), because the first term will be a large number of the order ${\displaystyle \frac{1}{13}}\times {10}^3$, while the second term is a very small number of the order $1.4\times {10}^{-8}$. Other values for the parameter $\beta$ have been obtained in the literature. For example, in the review article by J. Pascual Sanchez [23], it was obtained for the approximate radius of the satellite orbit $r_s=26,561$ [km] that

(47)${\displaystyle \frac{G_{\oplus }{M}_{\oplus }}{r_s{c}^2}}=0.167.{10}^{-9}.$

If multiplied by 2 (in order to obtain the constant $\beta =$${\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^{2}a}}$), this will give $0.334\times {10}^{-9}$, which is in full accord (up to the second digit) with the estimate in this paper $\beta =0.33\times {10}^{-9}$.

4. Definitions of Elliptic Integrals of Higher Order and Comparisonwith Some Statements from Standard Textbooks

4.1. Higher-Order Elliptic Integrals—Basic Definitions

It is necessary for this definition to be given because, further, when calculating the second part of the propagation time, we shall encounter elliptic integrals of the fourth order.

According to the general definition, elliptic integrals are of the type

(48)$\int R(y,\sqrt{G\left(y\right)})dy\mathrm{or}\int R(x,\sqrt{P\left(x\right)})dx,$

where $R(y,\tilde{y})$ and $R(x,\tilde{x})$ are rational functions of the variables $y,$$\tilde{y}$ or $x,$$\tilde{x}$ and $\tilde{y}=\sqrt{G\left(y\right)}$, $\tilde{x}=\sqrt{P\left(x\right)}$ are square roots of the arbitrary polynomials $G\left(y\right)$ and $P\left(x\right)$ of the fourth or of the third degree, respectively. In accord with the notations in the monograph [34], the fourth-order and the third-order polynomials will be of the form

${\tilde{y}}^2=G\left(y\right)=a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4,$

(49)${\tilde{x}}^2=P\left(x\right)=a{x}^3+b{x}^2+cx+d.$

So, in the general case, elliptic integrals of the n-th order and of the first kind are defined as

$J_n^{\left(4\right)}\left(y\right)=\int {\displaystyle \frac{y^{n}dy}{\sqrt{a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4}}},$

(50)$J_n^{\left(3\right)}\left(x\right)=\int {\displaystyle \frac{x^{n}dx}{\sqrt{a{x}^3+b{x}^2+cx+d}}}.$

Elliptic integrals of the n-th order and of the third kind are defined as

$H_n^{\left(4\right)}\left(y\right)=\int {\displaystyle \frac{dy}{{(y-c_1)}^n\sqrt{a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4}}},$

(51)$H_n^{\left(3\right)}\left(x\right)=\int {\displaystyle \frac{dx}{{(x-c_2)}^n\sqrt{a{x}^3+b{x}^2+cx+d}}}.$

In the above definitions, the upper indices “$\left(3\right)$” or “$\left(4\right)$” denote the order of the algebraic polynomial under the square root in the denominator, while the lower indice “n” denotes the order of the elliptic integral, related to the “$x^n$” or “$y^n$” terms in the nominator of (50) and in the denominator of (51) (for n-th order elliptic integrals of the third kind).

4.2. Comparison with Some Definitions from Standard Textbooks

In the textbook [41], a statement is expressed (although without any proof) that “higher-order elliptic integrals in the Legendre form can be expressed by means of the three standard integrals

$\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}},\int {\displaystyle \frac{y^{2}dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}},$

(52)$\int {\displaystyle \frac{dy}{(1+h{y}^2)\sqrt{(1-y^2)(1-q^2{y}^2)}}},$

where the parameter h can be also a complex number”. As asserted in [41], these integrals “cannot be expressed through elementary functions”. This assertion is not precise, because in the preceding sections, an example was found, when the second integral in (52) in fact can be expressed through an elementary function—for this case, by the logarithmic term in Equation (73). A similar statement to the one in [41] is given also in the monograph [38]: “Integrals of the type (48) $\int R(y,\sqrt{G\left(y\right)})dy$ and $\int R(x,\sqrt{P\left(x\right)})dx$ cannot be expressed by elementary functions in their final form, even if an extended understanding of this notion is considered”.

In the contemporary monograph [34], a statement is expressed in a correct and precise manner, namely, “Every elliptic integral can be represented in the form of a linear combination of a rational function of the variables y and $\tilde{y}$, integrals of a rational function of the variable y and also the integrals

$\int {\displaystyle \frac{dy}{\sqrt{G\left(y\right)}}},\int {\displaystyle \frac{ydy}{\sqrt{G\left(y\right)}}},$

(53)$\int {\displaystyle \frac{y^{2}dy}{\sqrt{G\left(y\right)}}},\int {\displaystyle \frac{dy}{(y-c)\sqrt{G\left(y\right)}}}”.$

In this statement, it is not mentioned that elliptic integrals cannot be expressed by elementary functions. One new analytical algorithm for expressing elliptic integrals in the Legendre form will be proposed in the next sections. Since the monographs [38,41] are older monographs, this example clearly shows that some basic facts about elliptic functions and integrals cannot be considered to be established forever and the theory may change in time. One basic and contemporary monograph on elliptic functions is [42] (see, especially, ch. 8, dedicated to elliptic integrals). A short summary of the various types of elliptic integrals and functions is given in Appendix A of the book [43], as well as applications of elliptic functions in some problems of biophysics.

4.3. Definitions of Elliptic Integrals of the First, of the Secondand of the Third Kind and of Zero Order

These integrals will be defined, because they were already encountered in the calculations of the propagation time of the signal, emitted by a satellite moving along a plane elliptical orbit.

The elliptic integrals of the first, second and third kind and of the zero order (see ch. 8 of the monograph by Armitage and Eberlein [42] and also [34]) are partial cases of the more generally defined integrals (50) and (51) The three types of integrals are

(54)$\int {\displaystyle \frac{dx}{\sqrt[.]{\left(1-x^2\right)\left(1-q^2{x}^2\right)}}}=\underset{0}{\overset{\varphi }{\int }}{\displaystyle \frac{d\varphi }{\sqrt[.]{1-q^2{sin}^2\varphi }}}=F(q,\varphi ),$

(55)$\int {\displaystyle \frac{\sqrt[.]{1-q^2{x}^2}}{\sqrt[.]{1-x^2}}}dx=\int \sqrt[.]{1-q^2{sin}^2\varphi }d\varphi =E(q,\varphi ),$

(56)$\int {\displaystyle \frac{dx}{(1+n{x}^2)\sqrt[.]{\left(1-x^2\right)\left(1-q^2{x}^2\right)}}}=\int {\displaystyle \frac{d\varphi }{(1+n{sin}^2\varphi )\sqrt[.]{1-q^2{sin}^2\varphi }}}=\Pi (q,\varphi ).$

In (54), $F(q,\varphi )$ is the commonly accepted notation for an elliptic integral of the first kind, in (55) $E(q,\varphi )$ is the notation for an elliptic integral of the second kind and, in (56), $\Pi (q,\varphi )$ is the notation for the elliptic integral of the third kind. In the above equalities, the transformation $x=sin\varphi$ has been used.

The contemporary definition of the elliptic integrals (50) for $J_n^{\left(4\right)}\left(y\right)$ and $J_n^{\left(3\right)}\left(x\right)$ and (51) for $H_n^{\left(4\right)}\left(y\right)$ and $H_n^{\left(3\right)}\left(x\right)$ in the monograph [34] is fully consistent with the definition in older monographs, such as [44].

4.4. Can Elliptic Integrals Be Considered to beAbelian Ones?

The elliptic integrals (50) and (51) cannot be considered to be a generalization of the s.c. abelian integral,

because after setting up

(57)$y={\displaystyle \frac{\sqrt[.]{1-q^2{x}^2}}{\sqrt[.]{1-x^2}}},$

(58)$y=(1+n{x}^2)\sqrt[.]{\left(1-x^2\right)\left(1-q^2{x}^2\right)},$

it can easily be seen that these algebraic curves are different from the curves (31) and (32). Consequently, the elliptic integrals (54)–(56) are not abelian integrals. Similarly, it can be proved that the found integral (29)

$T={\displaystyle \frac{a}{c}}\int \sqrt[.]{{\displaystyle \frac{a_1{y}^3+a_2{y}^2+a_{3}y}{b_1{y}^3+b_2{y}^2+b_{3}y+b_4}}}dy$

is not also an abelian integral. Again, it should reminded that the above integral is found as a consequence of the approximation $\beta ={\displaystyle \frac{2V}{c^2}}\ll 1$, which physically means that the gravitational potential is small.

5. Mathematical Structure of the Expression for the Propagation Timeof the Signal, Emitted by a Satellite on a Plane Elliptical Orbit: ThePropagation Time as a Sum of Zero-Order Elliptic Integrals of the First,Second and Third Rank

The first term in expression (41), after performing the simple transformation ${\displaystyle \frac{\pi }{2}}-E=\overline{E}$, can be brought to the following form:

$T_1=\underset{0}{\overset{E}{\int }}\sqrt[.]{1-e^2{cos}^{2}E}dE=\underset{0}{\overset{E}{\int }}\sqrt[.]{1-e^2{sin}^2({\displaystyle \frac{\pi }{2}}-E)}dE$

(59)$=-\underset{{\displaystyle \frac{\pi }{2}}}{\overset{{\displaystyle \frac{\pi }{2}}-\overline{E}}{\int }}\sqrt[.]{1-e^2{sin}^2\overline{E}}d\overline{E}.$

This represents an elliptic integral of the second kind. Note that the s.c.modulus of the elliptic integral is equal to the eccentricity e of the orbit.

The second term in Equation (41) after performing the substitution

(60)$\sqrt[.]{{\displaystyle \frac{1+ecosE}{1-ecosE}}}=\overline{y}$

and introducing the notations

(61)${\tilde{k}}^2={\displaystyle \frac{1-e}{1+e}}=q,{\displaystyle \frac{\overline{y}}{\tilde{k}}}=y,{\left({\displaystyle \frac{\overline{y}}{\tilde{k}}}\right)}^2={\tilde{y}}_1,{\tilde{y}}_1=-\tilde{\tilde{y}}$

can be written as a sum of two integrals

(62)$T_2=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}\int \sqrt[.]{{\displaystyle \frac{1+ecosE}{1-ecosE}}}dE=I_2^{A)}+I_2^{\left(B\right)}.$

It should be noted that $\tilde{k}$ in Equation (61) is simply a notation, representing the real expression $\tilde{k}=\sqrt{{\displaystyle \frac{1-e}{1+e}}}=\sqrt{q}$, since $0<{\displaystyle \frac{1-e}{1+e}}<1$ (the eccentricity of the orbit e is always less than one, $0<e<1$). We shall call this second part (62) of the expression for the propagation time the $O\left(c^{-3}\right)$ correction.

5.1. Mathematical Proof of the Real-Valuedness of the SeparateExpressions for the Propagation Time for the Case of Plane Elliptic Orbits

5.1.1. The First Real-Valued Elliptic Integral $I_2^{\left(A\right)}$ of theFirst Kind in the Weierstrass Form in the $O\left(c^{-3}\right)$ Correction of thePropagation Time

The first term $I_2^{\left(A\right)}$ in Equation (62) can be represented as

(63)$I_2^{\left(A\right)}={\displaystyle \frac{4GM}{c^3}}{\displaystyle \frac{1}{\tilde{k}\sqrt[.]{e^2-1}}}\int {\displaystyle \frac{d{\tilde{y}}_1}{\sqrt[.]{{\tilde{y}}_1\left({\tilde{y}}_1-\frac{1}{{\tilde{k}}^4}\right)\left({\tilde{y}}_1-1\right)}}}.$

In terms of the variable

(64)${\tilde{y}}_1=\left({\displaystyle \frac{1+e}{1-e}}\right).\left({\displaystyle \frac{1+ecosE}{1-ecosE}}\right),$

the integral (63) represents an elliptic integral of zero order and of the first kind, written in the Weierstrass form (i.e., a polynomial of third order under the square root in the denominator). Due to the presence of the $\sqrt[.]{e^2-1}$ term in the denominator (again, $0<e\le 1$), it might seem that expression (63) is imaginary. However, this is not true, because it can be rewritten as

(65)$I_2^{\left(A\right)}={\displaystyle \frac{4GM}{c^3}}{\displaystyle \frac{1}{\tilde{k}i\sqrt[.]{1-e^2.}\left(-i\right)}}\int {\displaystyle \frac{d\tilde{\tilde{y}}}{\sqrt[.]{\tilde{\tilde{y}}\left(\tilde{\tilde{y}}+1\right)\left(\tilde{\tilde{y}}+\frac{1}{{\tilde{k}}^4}\right)}}},$

thus representing an elliptic integral of zero order and of the first kind. The classification of elliptic integrals, depending on the degree of the polynomial in the denominator, is given in the monographs [34,38,45]. Since $\tilde{k}$ is a real-valued expression, the two imaginary units i in the denominator appear from $\sqrt[.]{e^2-1}=i\sqrt[.]{1-e^2}$ and also from the square root, after performing the variable transformation from ${\tilde{y}}_1$ to $\tilde{\tilde{y}}$ in ${\tilde{y}}_1\left({\tilde{y}}_1-{\displaystyle \frac{1}{{\tilde{k}}^4}}\right)\left({\tilde{y}}_1-1\right)$ in Equation (63). As a result, a factor $\sqrt{{(-1)}^3}=\sqrt{i^6}=i^3=-i$ will appear in the denominator of Equation (65). The whole expression Equation (65) turns out to be a real-valued one and with a negative sign, because $d{\tilde{y}}_1=-d\tilde{\tilde{y}}$. Since this expression is a part of the expression for the propagation time, this is physically reasonable.

5.1.2. The Second Real-Valued Elliptic Integral $I_2^{\left(B\right)}$ of theThird Kind

The second term $I_2^{\left(B\right)}$ in Equation (62) can also be written in the form of a real-valued expression, as it should be

(66)$I_2^{\left(B\right)}={\displaystyle \frac{4GM}{c^3{q}^2}}{\displaystyle \frac{1}{\sqrt[.]{1-e^2}}}\int {\displaystyle \frac{d{\tilde{y}}_1}{({\tilde{y}}_1-\frac{1}{q})\sqrt[.]{{\tilde{y}}_1({\tilde{y}}_1+1)({\tilde{y}}_1+\frac{1}{q^2})}}}.$

This is an elliptic integral of the third kind. Consequently, the whole expression for the propagation time can be represented as a sum of elliptic integrals of the second, first and third kind [34].

6. Propagation Time for a Signal, Emitted and Intercepted bySatellites on a Space-Oriented Orbit

6.1. Three-Dimensional Orbit Parametrization and the General Formulaefor the Orbit Parametrization

The formulaes for the three-dimensional parametrization of space-oriented orbits previously were given by the expressions (9)–(11). Again, after using the null cone Equation (40) $T=\int {\displaystyle \frac{v}{c}}\sqrt[.]{{\displaystyle \frac{(1-\frac{2V}{c^2})}{(1+\frac{2V}{c^2})}}}$$dt\approx \int {\displaystyle \frac{v}{c}}(1-{\displaystyle \frac{2V}{c^2}})dt$, remembering relation (15) $vd{t}_{cel}=v_{f}df$ and also expression (16) $v_f={\displaystyle \frac{na}{\sqrt{1-e^2}}}\sqrt{1+e^2+2ecosf}$ for the velocity $v_f$, associated with the true anomaly angle f, the following expression can be obtained for the propagation time $\overline{T}$ for the three-dimensional case:

(67)$\overline{T}=\int {\displaystyle \frac{v}{c}}(1-{\displaystyle \frac{2V}{c^2}})dt={\tilde{T}}_1+{\tilde{T}}_2={\displaystyle \frac{1}{c}}\int vdt-{\displaystyle \frac{2}{c^3}}\int vVdt.$

We denote the propagation time for the three-dimensional case (9)–(11) by $\overline{T}$, in order to distinguish it from the propagation time T for the case of a satellite, moving along a plane elliptical orbit.

6.2. Analytical Calculation of the First Elliptic Integral (First $O\left({\displaystyle \frac{1}{c}}\right)$ Correction) Without the Use of Elliptic Integrals

Let us take the first $O\left({\displaystyle \frac{1}{c}}\right)$ propagation time correction

(68)${\tilde{T}}_1={\displaystyle \frac{1}{c}}\int v_{f}df={\displaystyle \frac{na}{c\sqrt{1-e^2}}}\int \sqrt{1+e^2+2ecosf}df={\displaystyle \frac{na}{c\sqrt{1-e^2}}}{\tilde{T}}_1^{\left(1\right)}$

and let us perform the series of six subsequent variable transformations

(69)$\sqrt{1+e^2+2ecosf}=y,\overline{y}={(e+1)}^2-y^2,$

(70)$\sqrt{1-{\displaystyle \frac{\overline{y}{q}^2}{\overline{y}+B_1}}}=z,z^2=z_1,1-z_1=m,\overline{m}=m-{\displaystyle \frac{q^2}{2}}.$

The first transformation $\sqrt{1+e^2+2ecosf}=y$ in (69) will be used only in this paragraph; further, another transformation $\tilde{y}(1+e)=\sqrt{1+e^2+2ecosf}$ will be applied.

After applying the sequence of transformations (69) and (70), the following sum of two integrals for ${\tilde{T}}_1$ will be obtained:

(71)${\tilde{T}}_1=-{\displaystyle \frac{nai}{c\sqrt{(1-e^2)q}}}\left[\int {\displaystyle \frac{d\overline{m}}{\left(\overline{m}-\frac{q^2}{2}\right)\sqrt{\overline{m}-\frac{q^2}{2}}}}+\int {\displaystyle \frac{d\overline{m}}{\left(\overline{m}+\frac{q^2}{2}\right)\sqrt{\overline{m}-\frac{q^2}{2}}}}\right],$

where, again, the notation $q={\displaystyle \frac{1-e}{1+e}}$ has been used and each of the integrals inside the square bracket will be denoted correspondingly as ${\tilde{T}}_1^{\left(1\right)}$ and ${\tilde{T}}_1^{\left(2\right)}$. Making use of the analytically calculated integral from the book [46],

$\int {\displaystyle \frac{dx}{(x\pm p)\sqrt{a+2bx+c{x}^2}}}$

(72)$=-{\displaystyle \frac{1}{\sqrt{a\mp 2bp+c{p}^2}}}ln\left\{{\displaystyle \frac{\sqrt{a+2bx+c{x}^2}+\sqrt{a\mp 2bp+c{p}^2}}{x\pm p}}+{\displaystyle \frac{b\mp cp}{\sqrt{a\mp 2bp+c{p}^2}}}\right\},$

one can derive for the second integral ${\tilde{T}}_1^{\left(2\right)}$ in Equation (71) the following expression:

(73)${\tilde{T}}_1^{\left(2\right)}=-{\displaystyle \frac{na(1+e)\sqrt{2}}{c\sqrt{1-e^2}\sqrt{q}\sqrt{3e^2+2e+3}}}ln\left[\left({\displaystyle \frac{\sqrt{2}{m}_1(f_b;r_b)+\frac{1}{2}{m}_2(f_b;r_b)}{\sqrt{2}{m}_1(f_a;r_a)+\frac{1}{2}{m}_2(f_a;r_a)}}\right)\left({\displaystyle \frac{m_3(f_a;r_a)}{m3(f_b;r_b)}}\right)\right].$

In (73), $m_1(f;r)$, $m_2(f;r)$, $m_3(f;r)$ are expressions, written in terms either of the initial and final true anomaly angles $f_a$ and $f_b$ or of the initial distance $r_a$ (at which the emission of the signal takes place) and the final point $r_b$ of reception, corresponding to the propagation time ${\tilde{T}}_1^{\left(2\right)}$. The first integral ${\tilde{T}}_1^{\left(1\right)}$ can be calculated analogously. It should be noted that both ${\tilde{T}}_1^{\left(1\right)}$ and ${\tilde{T}}_1^{\left(2\right)}$ are real-valued expressions and moreover, the logarithmic term is again present, as in the original Shapiro delay formulae (1).

Further, this analytical method will not be used because a transformation will be proposed, enabling to find analytical expressions for zero-order elliptic integrals in the Weierstrass or in the Legendre form.

6.3. Analytical Calculation of the First $O\left({\displaystyle \frac{1}{c}}\right)$ Correctionby Means of Elliptic Integrals

Let us calculate the integral (68) ${\tilde{T}}_1={\displaystyle \frac{na}{c\sqrt{1-e^2}}}\int \sqrt{1+e^2+2ecosf}df$ by making use of the substitution

(74)$y=\sqrt{{\displaystyle \frac{1}{q}}{\displaystyle \frac{(1+ecosE)}{(1-ecosE)}}},q={\displaystyle \frac{1-e}{1+e}}$

and also of the well-known relation from celestial mechanics between the eccentric anomaly angle E and the true anomaly angle f [18]:

(75)$tan{\displaystyle \frac{f}{2}}=\sqrt{{\displaystyle \frac{1-cosf}{1+cosf}}}=\sqrt{{\displaystyle \frac{1+e}{1-e}}}tan{\displaystyle \frac{E}{2}}.$

Then, we can write the integral ${\tilde{T}}_1$ in the form of an elliptic integral of the second order and of the first kind in the Legendre form

(76)${\tilde{T}}_1=-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}\int {\displaystyle \frac{y^{2}dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}=-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}{\tilde{J}}_2^{\left(4\right)}(y;q).$

Thus, the integral (68) is very interesting from a mathematical point of view, because it becomes evident that it is an elliptic one after performing the transformations (74) and (75). The integral ${\tilde{J}}_2^{\left(4\right)}(y;q)$ in (76) can also be represented as

(77)${\tilde{J}}_2^{\left(4\right)}(y;q)=\int {\displaystyle \frac{y^{2}dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}$

(78)$=-{\displaystyle \frac{1}{q^2}}\int {\displaystyle \frac{(1-q^2{y}^2)dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}+{\displaystyle \frac{1}{q^2}}\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}$

(79)$=-{\displaystyle \frac{1}{q^2}}\int \sqrt{1-q^2{sin}^2\phi }d\phi +{\displaystyle \frac{1}{q^2}}\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}.$

The first integral in Equation (79) is an elliptic integral of the second kind (denoted usually by $E\left(\phi \right)=\int \sqrt{1-q^2{sin}^2\phi }d\phi$, where $x=sin\phi$). So we obtain a relation between the second-order elliptic integral ${\tilde{J}}_2^{\left(4\right)}(y;q)$ of the first kind in the Legendre form, the zero-order elliptic integral of the first kind in the Legendre form ${\tilde{J}}_0^{\left(4\right)}(y;q)$ (the second integral in (79)) and the elliptic integral $E\left(\phi \right)$ (see the monograph [34]):

(80)${\tilde{J}}_2^{\left(4\right)}(y;q)=\int {\displaystyle \frac{y^{2}dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}=-{\displaystyle \frac{1}{q^2}}E\left(\phi \right)+{\displaystyle \frac{1}{q^2}}{\tilde{J}}_0^{\left(4\right)}(y;q).$

Since, in the preceding section, an analytical expression was obtained for ${\tilde{T}}_1$ without the use of any elliptic functions, the obtained result in Equation (79) means that for the investigated case, elliptic integrals of second order can be expressed through elementary functions.

6.4. Second Analytical Calculation of the First Time $O\left({\displaystyle \frac{1}{c}}\right)$Correction in Terms of Second-Order Elliptic Integrals and Proof of theReal-Valuedness of the $O\left({\displaystyle \frac{1}{c}}\right)$ Correction

This second calculation will not make any use of the eccentric anomaly angle variable E. The main purpose of this second analytical calculation will be prove the real-valuedness of the expression (76)

${\tilde{T}}_1=-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}\int {\displaystyle \frac{y^{2}dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}=-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}{J}_2^{\left(4\right)}(y,q),$

which means that the second-order elliptic integral in the Legendre form should be imaginary-valued.

Let us apply the variable transformation

(81)$\tilde{y}={\displaystyle \frac{\sqrt{1+2ecosf+e^2}}{1+e}}$

in expression (68) for ${\tilde{T}}_1$, after which the integral ${\tilde{T}}_1$ acquires the form

(82)${\tilde{T}}_1={\displaystyle \frac{na}{c\sqrt{1-e^2}}}\int \sqrt{1+e^2+2ecosf}df$

$=i{\displaystyle \frac{2na(1+e)}{cq\sqrt{1-e^2}}}\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-\frac{{\tilde{y}}^2}{q^2}\right)}}}$

(83)$={\displaystyle \frac{2na(1+e)}{cq\sqrt{1-e^2}}}\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(\frac{{\tilde{y}}^2}{q^2}-1\right)}}}$

(84)$={\displaystyle \frac{2na(1+e)}{icq\sqrt{1-e^2}}}\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-\frac{{\tilde{y}}^2}{q^2}\right)}}}={\displaystyle \frac{2na(1+e)}{icq\sqrt{1-e^2}}}{\tilde{J}}_2^{\left(4\right)}(\tilde{y};{\displaystyle \frac{1}{q}})$

Note that the coefficient in front of the integral is modified in comparison with the one in Equation (76). More importantly, the resulting expression is again a real-valued one due to the property of the elliptic integral

(85)${\tilde{J}}_2^{\left(4\right)}(\tilde{y},{\displaystyle \frac{1}{q}})=-i\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(\frac{{\tilde{y}}^2}{q^2}-1\right)}}},$

which is easily established after comparing expressions (83) and (84). So, it can be concluded that ${\tilde{T}}_1$, in fact, is expressed in two different variables. The first one is (76)

${\tilde{T}}_1=-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}\int {\displaystyle \frac{y^{2}dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}=-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}{\tilde{J}}_2^{\left(4\right)}(y;q),$

where the variable y is expressed as (74) $y=\sqrt{{\displaystyle \frac{1}{q}}{\displaystyle \frac{(1+ecosE)}{(1-ecosE)}}}$, $q={\displaystyle \frac{1-e}{1+e}}$, and the second one is (83) in terms of the variable (81) $\tilde{y}={\displaystyle \frac{\sqrt{1+2ecosf+e^2}}{1+e}}$.

Formulae (85) also is valid, because due to the inequality $cosf\le 1$ and the choice of the variable $\tilde{y}$ in Equation (81), it can easily be proved that

(86)$\tilde{y}\le 1⟹1-{\tilde{y}}^2\ge 0,1-{\displaystyle \frac{{\tilde{y}}^2}{q^2}}\ge -{\displaystyle \frac{(1-q^2)}{q^2}}=-{\displaystyle \frac{4e}{{(1-e)}^2}}.$

Consequently, since $1-{\displaystyle \frac{{\tilde{y}}^2}{q^2}}$ can take negative values (note that $q^2={\left({\displaystyle \frac{1-e}{1+e}}\right)}^2<1$) and, thus, $0<{\displaystyle \frac{{\tilde{y}}^2}{q^2}}-1\le {\displaystyle \frac{(1-q^2)}{q^2}}$, the representation (85) gives a real-valued integral $\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(\frac{{\tilde{y}}^2}{q^2}-1\right)}}}$ and, thus, an imaginary integral ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},{\displaystyle \frac{1}{q}})$ is obtained.

Let us present another more elegant and simple proof (see also the review paper [22]) that the inequalities (86) are fulfilled. Inverting the sign of the last inequality in (86), the inequality

(87)${\displaystyle \frac{{\tilde{y}}^2}{q^2}}-1\le {\displaystyle \frac{4e}{{(1-e)}^2}}$

should be proved. Keeping in mind the definition for the variable $\tilde{y}$ (81) $\tilde{y}={\displaystyle \frac{\sqrt{1+2ecosf+e^2}}{1+e}}$ and also for $q^2$; the inequality (87) is transformed to

${\displaystyle \frac{(1+2ecosf+e^2)}{{(1+e)}^2}}-{\displaystyle \frac{{(1-e)}^2}{{(1+e)}^2}}$

(88)$\le {\displaystyle \frac{4e}{{(1-e)}^2}}{q}^2={\displaystyle \frac{4e}{{(1+e)}^2}}.$

From here, it follows that

(89)$1+2ecosf+e^2\le 4e+{(1-e)}^2\Rightarrow 2ecosf\le 2e,$

which is fulfilled for all eccentricities e, because of the simple trigonometric inequality $cosf\le 1$ for the function cos and also because $0<e<1$. Therefore, inequality (87) is fulfilled.

The other case, which has to be considered, is when $1-{\displaystyle \frac{{\tilde{y}}^2}{q^2}}$ takes only positive values:

(90)$1-{\displaystyle \frac{{\tilde{y}}^2}{q^2}}\ge 0\ge -{\displaystyle \frac{(1-q^2)}{q^2}}$

and, consequently, ${\displaystyle \frac{{\tilde{y}}^2}{q^2}}-1\le 0$. However, in view of the future calculations, we will be interested in the case when the transformation $q\to {\displaystyle \frac{1}{q}}$ is applied in (85), so that the integral ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$ is obtained. Then, from (90) after the transformation $q\to {\displaystyle \frac{1}{q}}$, it follows that ${\tilde{y}}^2{q}^2-1\le 0$ or $1-{\tilde{y}}^2{q}^2\ge 0$. Thus, one will have the integral transformation

(91)$\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(\frac{{\tilde{y}}^2}{q^2}-1\right)}}}\stackrel{q\to {\displaystyle \frac{1}{q}}}{⟹}\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}}},$

where the integral in the right is real-valued. Multiplying both integrals by ${\displaystyle \frac{1}{i}}$ and remembering the definition (85) for ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},{\displaystyle \frac{1}{q}})$, one obtains

(92)${\tilde{J}}_2^{\left(4\right)}(\tilde{y},{\displaystyle \frac{1}{q}})\stackrel{q\to {\displaystyle \frac{1}{q}}}{⟹}{\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)={\displaystyle \frac{1}{i}}\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}}}.$

This means that ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$ is imaginary valued. After presenting the calculation for the second $O\left(c^{-3}\right)$ propagation time correction in the next sections, it will become clear that due to the important transformation (92), it will be possible to prove the real-valuedness of terms with second-order elliptic integrals in the final expression for the propagation time for the case of space-oriented orbits.

6.5. An Earlier Example for an Elliptic Integral, Expressed AlsoAnalytically

However, from the formulaes in the preceding section, it cannot be concluded that only abelian integrals can be treated analytically (in the sense of being expressed by elementary functions) and the elliptic integrals cannot be expressed by analytical functions. A typical example was given earlier, related to the calculation of the integral (68)

${\tilde{T}}_1={\displaystyle \frac{na}{c\sqrt{1-e^2}}}\int \sqrt{1+e^2+2ecosf}df,$

which was proved also to be the elliptic integral (83) ${\tilde{T}}_1={\displaystyle \frac{2na(1+e)}{cq\sqrt{1-e^2}}}\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(\frac{{\tilde{y}}^2}{q^2}-1\right)}}}$ of second order in the Legendre form in terms of the variable (81) $\tilde{y}={\displaystyle \frac{\sqrt{1+2ecosf+e^2}}{1+e}}$.

Also, in terms of the variables (74)

$y=\sqrt[.]{{\displaystyle \frac{(1+ecosE)}{q(1-ecosE)}}},(q={\displaystyle \frac{1-e}{1+e}})$

and (75)

$tan{\displaystyle \frac{f}{2}}=\sqrt[.]{{\displaystyle \frac{1-cosf}{1+cosf}}}=\sqrt[.]{{\displaystyle \frac{1+e}{1-e}}}tan{\displaystyle \frac{E}{2}},$

the integral ${\tilde{T}}_1$ was expressed as an elliptic integral of the second order and in the Legendre form (76)

${\tilde{T}}_1=-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}\int {\displaystyle \frac{y^{2}dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}=-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}{\tilde{J}}_2^{\left(4\right)}(y;q).$

Earlier, by means of the transformations (69) and (70), it was proved that the same integral can be expressed analytically by means of the integrals (71)

${\tilde{T}}_1=-{\displaystyle \frac{nai}{c\sqrt{(1-e^2)q}}}\left[\int {\displaystyle \frac{d\overline{m}}{\left(\overline{m}-\frac{q^2}{2}\right)\sqrt{\overline{m}-\frac{q^2}{2}}}}+\int {\displaystyle \frac{d\overline{m}}{\left(\overline{m}+\frac{q^2}{2}\right)\sqrt{\overline{m}-\frac{q^2}{2}}}}\right],$

which can be integrated. Thus, it can be concluded that when expressed in two different variables $\tilde{y}$ and $\overline{m}$, one and the same integral (68) can be considered to be an elliptic one, but also, it can be written in another form and even calculated exactly. From the transformations (69) and (70), the relation between the variables f and $\overline{m}$ can be found to be

(93)$cosf=1+{\displaystyle \frac{1}{2e}}{\displaystyle \frac{(2\overline{m}+q^2)B_1}{(2\overline{m}-q^2)}}.$

If account is taken also of (81) for $\tilde{y}$, the relation between the variables $\tilde{y}$ and $\overline{m}$ is calculated to be a quadratic algebraic curve with respect to the variable $\tilde{y}$:

(94)${\tilde{y}}^2=1+{\displaystyle \frac{1}{{(1+e)}^2}}.{\displaystyle \frac{(2\overline{m}+q^2)B_1}{(2\overline{m}-q^2)}}.$

Then, if the integral (68) in terms of the variable $\overline{m}$ can be integrated, it can be asked if there is any motivation for the assertion in the monographs [38,41,42] that elliptic integrals cannot be integrated in elementary functions. Further, it shall be proved that analytical formulaes can be found for the determination of the s.c. Weierstrass invariants $g_2$ and $g_3$, depending on complicated polynomials of the modulus q of the elliptic integral in the Legendre form. But, then, the s.c. “four-dimensional parametrization” of an elliptic curve will allow the parametrization of an fourth-order polynomial in terms of expressions, depending ot the Weierstrass function $\rho \left(z\right)$ and its derivative ${\displaystyle \frac{d\rho \left(z\right)}{dz}}$, which, of course, are not elementary functions. This will be proved further, but evidently, depending on the variables used in the calculations, a certain elliptic integral can be expressed analytically in elementary functions, but also can be expressed in an analytical form, containing non-elementary functions.

7. Fourth-Order Elliptic Integrals and the Second $\mathit{O}\left({\displaystyle \frac{\mathbf{1}}{{\mathit{c}}^{\mathbf{3}}}}\right)$Correction of the Propagation Time for the Case of a Space-Oriented Orbit

7.1. The First Part of the Second $O\left({\displaystyle \frac{1}{c^3}}\right)$ Correction,Expressed in Terms of Elliptic Integrals of the Second Order

Taking into account the general equality (67) $\overline{T}={\tilde{T}}_1+{\tilde{T}}_2={\displaystyle \frac{1}{c}}\int vdt-{\displaystyle \frac{2}{c^3}}\int vVdt$, the second $O\left({\displaystyle \frac{1}{c^3}}\right)$ correction can be written as

(95)${\tilde{T}}_2=-\int {\displaystyle \frac{v}{c}}.{\displaystyle \frac{2V}{c^2}}dt=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}\int {\displaystyle \frac{v_f}{r}}df,$

where the radius-vector r in the orbital plane, the velocity $v_f$, associated with the true anomaly angle f and the potential V of the gravitational field around the Earth are defined in the usual way:

$v_f={\displaystyle \frac{na}{\sqrt{1-e^2}}}\sqrt{1+2ecosf+e^2},r={\displaystyle \frac{a(1-e^2)}{1+ecosf}},$

(96)$V={\displaystyle \frac{GM}{r}}={\displaystyle \frac{GM}{a(1-e^2)}}(1+ecosf).$

The final integral for the second correction is of the form

(97)${\tilde{T}}_2=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}.{\displaystyle \frac{na}{a\sqrt{1-e^2}}}\int (1+ecosf)\sqrt{1+2ecosf+e^2}df=T_2^{\left(1\right)}+T_2^{\left(2\right)},$

where

$T_2^{\left(1\right)}=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}.{\displaystyle \frac{n}{{(1-e^2)}^{\frac{3}{2}}}}\int \sqrt{1+2ecosf+e^2}df$

(98)$=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}.{\displaystyle \frac{n}{{(1-e^2)}^{\frac{3}{2}}}}{\tilde{T}}_1.$

In (98), ${\tilde{T}}_1$ is the previously calculated integral ${\tilde{T}}_1=\int \sqrt{1+2ecosf+e^2}df$. The second term $T_2^{\left(2\right)}$ in the second correction ${\tilde{T}}_2$ in (97) is

$T_2^{\left(2\right)}=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}.{\displaystyle \frac{ne}{{(1-e^2)}^{\frac{3}{2}}}}\int cosf\sqrt{1+2ecosf+e^2}df$

(99)$=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}{\displaystyle \frac{ne}{{(1-e^2)}^{\frac{3}{2}}}}{\tilde{T}}_2^{\left(2\right)},$

where ${\tilde{T}}_2^{\left(2\right)}$ is the notation for the more complicated integral

(100)${\tilde{T}}_2^{\left(2\right)}=\int cosf\sqrt{1+2ecosf+e^2}df.$

7.2. The Second Part of the Second $O\left({\displaystyle \frac{1}{c^3}}\right)$ Correction,Expressed in Terms of Elliptic Integrals of the Second and of the FourthOrder

It can be proved that

(101)$T_2^{\left(2\right)}=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}.{\displaystyle \frac{ne}{{(1-e^2)}^{\frac{3}{2}}}}{\tilde{T}}_2^{\left(2\right)}$

$=-in{q}^{{\displaystyle \frac{5}{2}}}{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}(1+e^2){\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$

(102)$+in{q}^{{\displaystyle \frac{3}{2}}}{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}{\displaystyle \frac{(1+e^2)}{(1-e^2)}}{\tilde{J}}_4^{\left(4\right)}(\tilde{y},q),$

where

${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)=\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}}}$

(103)$={\displaystyle \frac{1}{i}}\int {\displaystyle \frac{q^3{\widehat{y}}^{2}d\widehat{y}}{\sqrt{\left({\widehat{y}}^2-1\right)\left(1-q^2{\widehat{y}}^2\right)}}}={\displaystyle \frac{q^3}{i}}{\tilde{J}}_2^{\left(4\right)}(\widehat{y},q)$

is an elliptic integral of the second order in the Legendre form and $\widehat{y}$ is the notation for the variable $\widehat{y}={\displaystyle \frac{\tilde{y}}{q}}$. The last formulae can also be compared to (85). Analogously, the formulae for the fourth-order elliptic integral ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$ can be written as

(104)${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)={\displaystyle \frac{q^5}{i}}\int {\displaystyle \frac{{\widehat{y}}^{4}d\widehat{y}}{\sqrt{\left({\widehat{y}}^2-1\right)\left(1-q^2{\widehat{y}}^2\right)}}}.$

7.3. Relation Between the Fourth-Order and the Second-Order EllipticIntegrals in the Expression for $T_2^{\left(2\right)}$

From the recurrent chain (20), after setting up $n=0$ and after integration, it can easily be established that

$\int \sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}d\tilde{y}={\displaystyle \frac{2}{3}}{\tilde{J}}_0^{\left(4\right)}(\tilde{y},q)$

(105)$-{\displaystyle \frac{1}{3}}(1+q^2){\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)+{\displaystyle \frac{1}{3}}\left[\tilde{y}\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}\right]{\mid }_{\tilde{y}={\tilde{y}}_0}^{\tilde{y}={\tilde{y}}_1}.$

We have already proved that ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$ can be expressed in elementary functions, yet we do not know whether ${\tilde{J}}_0^{\left(4\right)}(\tilde{y},q)$ or the second-rank integral

(106)$\int \sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}d\tilde{y}$

can also be expressed through elementary functions. This will be proved further.

Calculating the derivatives

(107)${\displaystyle \frac{d}{d\tilde{y}}}\left(\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}\right)\mathrm{and}{\displaystyle \frac{d}{d\tilde{y}}}\left(\tilde{y}\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}\right),$

integrating the resulting equations from ${\tilde{y}}_0$ to ${\tilde{y}}_1$, and combining all the three equations, the following two equations can be derived for ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$ and ${\tilde{J}}_3^{\left(4\right)}(\tilde{y},q)$:

${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)={\displaystyle \frac{1}{3q^2}}\left[\tilde{y}\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}\right]{\mid }_{\tilde{y}={\tilde{y}}_0}^{\tilde{y}={\tilde{y}}_1}$

(108)$+{\displaystyle \frac{2(1+q^2)}{3q^2}}{\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)-{\displaystyle \frac{1}{3q^2}}{\tilde{J}}_0^{\left(4\right)}(\tilde{y},q),$

(109)${\tilde{J}}_3^{\left(4\right)}(\tilde{y},q)={\displaystyle \frac{1}{2q^2}}\left[\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}\right]{\mid }_{\tilde{y}={\tilde{y}}_0}^{\tilde{y}={\tilde{y}}_1}+{\displaystyle \frac{(1+q^2)}{2q^2}}{\tilde{J}}_1^{\left(4\right)}.$

In expressions (105)–(109), the symbol “${\mid }_{\tilde{y}={\tilde{y}}_0}^{\tilde{y}={\tilde{y}}_1}$” means that the corresponding expression is taken at the value $\tilde{y}={\tilde{y}}_1$ of the upper boundary and, from it, the value of the expression at $\tilde{y}={\tilde{y}}_0$ (the value at the lower boundary) is subtracted.

We should also add to this system of equations the earlier derived equation for ${\tilde{J}}_2^{\left(4\right)}$:

${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)=\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{(1-{\tilde{y}}^2)(1-q^2{\tilde{y}}^2)}}}$

(110)$=-{\displaystyle \frac{1}{q^2}}E\left(\phi \right)+{\displaystyle \frac{1}{q^2}}{\tilde{J}}_0^{\left(4\right)}({\tilde{y}}^2;q).$

The second Equation (109) for ${\tilde{J}}_3^{\left(4\right)}$ clearly suggests that ${\tilde{J}}_3^{\left(4\right)}$ is expressed in elementary functions because the first-order elliptic integral ${\tilde{J}}_1^{\left(4\right)}$$=\int {\displaystyle \frac{\tilde{y}d\tilde{y}}{\sqrt{(1-{\tilde{y}}^2)(1-q^2{\tilde{y}}^2)}}}$ can be analytically calculated after performing the Euler substitution:

(111)$\sqrt{{\tilde{y}}^4-{\displaystyle \frac{(1+q^2)}{q^2}}{\tilde{y}}^2+{\displaystyle \frac{1}{q^2}}}=u+{\tilde{y}}^2.$

The result is

(112)${\tilde{J}}_1^{\left(4\right)}=-{\displaystyle \frac{1}{2q^3}}ln\mid \sqrt{(1-{\tilde{y}}^2)(1-q^2{\tilde{y}}^2)}-{\tilde{y}}^2+{\displaystyle \frac{(1+q^2)}{2q^2}}{\mid }_{\tilde{y}={\tilde{y}}_0}^{\tilde{y}={\tilde{y}}_1}.$

Concerning the fourth-order elliptic integral ${\tilde{J}}_4^{\left(4\right)}$, taking into account the already proved fact that ${\tilde{J}}_2^{\left(4\right)}$ can be expressed in elementary functions, it would follow that ${\tilde{J}}_4^{\left(4\right)}$ can also be expressed by elementary functions. However, it is necessary to prove that the zero-order elliptic integral ${\tilde{J}}_0^{\left(4\right)}$ can be expressed in elementary functions. Further, it will be proved that the integral can be expressed by an analytical formulae, but also, it can be numerically calculated, since the modulus of the elliptic integral is equal to the eccentricity $e=0.01323881349526$ of the GPS orbit, which is exactly known [13]. The general mathematical theory for the recurrent relations between n-th order elliptic integrals of any kind and the lower-order elliptic integrals has been developed in the monographs [34,38,45].

7.4. The Several Parts of the Entire Expression for the PropagationTime for the Case of Space-Oriented Orbits

Now, we shall collect all the separate terms, which constitute the expression for the propagation time T. These terms are (83), earlier proved to be a real-valued integral

${\tilde{T}}_1={\displaystyle \frac{2na(1+e)}{cq\sqrt{1-e^2}}}\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(\frac{{\tilde{y}}^2}{q^2}-1\right)}}},$

the expression (98) for $T_2^{\left(1\right)}$

$T_2^{\left(1\right)}=-{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}.{\displaystyle \frac{n}{{(1-e^2)}^{\frac{3}{2}}}}{\tilde{T}}_1,$

expression (102) for $T_2^{\left(2\right)}$

$T_2^{\left(2\right)}=-in{q}^{{\displaystyle \frac{5}{2}}}{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}(1+e^2){\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$

$+in{q}^{{\displaystyle \frac{3}{2}}}{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }}{c^3}}{\displaystyle \frac{(1+e^2)}{(1-e^2)}}{\tilde{J}}_4^{\left(4\right)}(\tilde{y},q),$

expression (103) for ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$

${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)=\int {\displaystyle \frac{{\tilde{y}}^{2}d\tilde{y}}{\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}}}$

$={\displaystyle \frac{q^3}{i}}\int {\displaystyle \frac{{\widehat{y}}^{2}d\widehat{y}}{\sqrt{\left({\widehat{y}}^2-1\right)\left(1-q^2{\widehat{y}}^2\right)}}}={\displaystyle \frac{q^3}{i}}{\tilde{J}}_2^{\left(4\right)}(\widehat{y},q),$

expression (104) for $J_4^{\left(4\right)}(\tilde{y},q)$

${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)={\displaystyle \frac{q^5}{i}}\int {\displaystyle \frac{{\widehat{y}}^{4}d\widehat{y}}{\sqrt{\left({\widehat{y}}^2-1\right)\left(1-q^2{\widehat{y}}^2\right)}}}.$

Note that this expression for ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$ resembles the previous expression ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$, since both of them have one and the same denominators $\sqrt{\left({\widehat{y}}^2-1\right)\left(1-q^2{\widehat{y}}^2\right)}$, which determine whether both formulaes are imaginary- or real- valued. Earlier, by formulaes (82)–(86) and (90)–(92), it was proved that ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$ is an imaginary-valued expression. Consequently, it can be concluded that expression (104) for ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$ is also imaginary-valued.

This conclusion can be made also from the chain of elliptic integrals, establishing the relation between elliptic integrals of zero, second and fourth order (108):

${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)={\displaystyle \frac{1}{3q^2}}\left[\tilde{y}\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}\right]{\mid }_{\tilde{y}={\tilde{y}}_0}^{\tilde{y}={\tilde{y}}_1}$

$+{\displaystyle \frac{2(1+q^2)}{3q^2}}{\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)-{\displaystyle \frac{1}{3q^2}}{\tilde{J}}_0^{\left(4\right)}(\tilde{y},q).$

The integrals ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$ and ${\tilde{J}}_0^{\left(4\right)}(\tilde{y},q)$ are imaginary-valued. This is related to the imaginary-valuedness of the expression $\sqrt{\left(1-{\tilde{y}}^2\right)\left(1-q^2{\tilde{y}}^2\right)}$ (the first term in (108)), which enters also the denominators of these two integrals. Then, the left-hand side of (108) (the integral ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$) will also be imaginary valued. This precludes the proof that ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$ is also imaginary-valued. Note also that if the variable $\tilde{y}$ in ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$ and ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$ is replaced by y, the integrals ${\tilde{J}}_2^{\left(4\right)}(y,q)$ and ${\tilde{J}}_4^{\left(4\right)}(y,q)$ will also be imaginary valued.

7.5. The Entire Expression for the Propagation Time of the Signal forthe Case of Space-Oriented Orbits and for All Orders in $1/c$

Making use of all the above mentioned formulaes, the total propagation time (both the $O\left({\displaystyle \frac{1}{c}}\right)$ and $O\left({\displaystyle \frac{1}{c^3}}\right)$ corrections) can be written as

$\overline{T}=\left[-2i{\displaystyle \frac{na}{c}}{q}^{{\displaystyle \frac{3}{2}}}+4i{\displaystyle \frac{G_{\oplus }{M}_{\oplus }{n}^{2}a}{c^4{(1-e^2)}^{\frac{3}{2}}}}{q}^{{\displaystyle \frac{3}{2}}}\right]{\tilde{J}}_2^{\left(4\right)}(y,q)$

(113)$-i{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }na{q}^{\frac{5}{2}}(1+e^2)}{c^3}}{\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)+i{\displaystyle \frac{2G_{\oplus }{M}_{\oplus }n{q}^{\frac{3}{2}}(1+e^2)}{c^3(1-e^2)}}{\tilde{J}}_4^{\left(4\right)}(\tilde{y},q),$

where all the integrals ${\tilde{J}}_2^{\left(4\right)}(y,q)$, ${\tilde{J}}_2^{\left(4\right)}(\tilde{y},q)$ and ${\tilde{J}}_4^{\left(4\right)}(\tilde{y},q)$ (74), written in terms of the variable $y=\sqrt{{\displaystyle \frac{1}{q}}{\displaystyle \frac{(1+ecosE)}{(1-ecosE)}}}$ or (81) $\tilde{y}={\displaystyle \frac{\sqrt{1+2ecosf+e^2}}{1+e}}$, are proved to be imaginary-valued.

This precludes the proof that the propagation time $\overline{T}$ is a real-valued expression, because all the imaginary-valued integrals are multiplied by imaginary coefficients.

8. Numerical Calculation of the Propagation Time of a Signal, Emittedby a Satellite on a Plane Elliptical Orbit

8.1. Numerical Calculation of the First Six Iterations of theEccentric Anomaly Angle E

We shall perform some numerical calculations of the propagation time, based on the derived formulae (41). The numerical data regarding the parameters of the $GPS$ orbit are taken from the PhD dissertation [13] and are known with great precision. From all the parameters, listed below, only the first three will be used in the calculation, performed by the online program web2.0 scientific calculator [47].

semi-major axis a	26,560.25169632944 [km]
eccentricity e	0.01323881349526,
mean anomaly M	−0.3134513508155 [rad],
inclination I	0.9614884100802 [rad],
longitude of the ascending node $\Omega$	−0.4495096737336 [rad],
argument of perigee $\omega$	−3.001488651204 [rad].

In the following numerical calculations, we shall use only the first three parameters $a,e$ and M, but in other publications, the approach will be extended and it will be required to find the propagation time dependent on the full set of six Kepler parameters $(a,e,M,I,\Omega ,\omega )$. Instead of M, either the eccentric anomaly angle E or the true anomaly angle f will be used. In fact, all the Kepler parameters might be used for finding the propagation time of the signal, and this is a problem to be solved in future.

For finding the propagation time for the case of a signal, emitted by a satellite, moving along a plane elliptical orbit with an eccentricity e, most important is to find the eccentric anomaly angle E from the Kepler equation $M=E-esinE$. The calculation will be based on formulaes $\left(2.54\right)$ and $\left(2.55\right)$ in Ch. 2 of the monograph [48].

The iterative sequence of the formulaes is

(114)$E_{(i+1)}=M+esin{E}_{\left(i\right)},i=0,1,2,$

where, for the first approximation, it is assumed $E_0=M=-0.3134513508155$ [rad]. The first three iterative solutions are given according to the following formulaes:

(115)$E_{\left(1\right)}=M+esinM,$

(116)$E_{\left(2\right)}=M+esin{E}_{\left(1\right)}=M+esin(M+esinM),$

(117)$E_{\left(3\right)}=M+esin{E}_{\left(2\right)}=M+esin[M+esin(M+esinM)].$

Thus, from the formulaes for $E_{\left(3\right)}$, $E_{\left(5\right)}$ and $E_{\left(6\right)}$, the following numerical values can be found:

(118)$E_{\left(3\right)}=M+esin{E}_{\left(2\right)}=-0.31758547588467897473\left[\mathrm{rad}\right],$

(119)$E_{\left(4\right)}=M+esin{E}_{\left(3\right)}=-0.31758548401096719083\left[\mathrm{rad}\right],$

(120)$E_{\left(5\right)}=M+esin{E}_{\left(4\right)}=-0.31758548411317083102\left[\mathrm{rad}\right],$

(121)$E_{\left(6\right)}=M+esin{E}_{\left(5\right)}=-0.31758548411445499592\left[\mathrm{rad}\right].$

The approximation $E_{\left(5\right)}$ up to the ninth digit is identical to $E_{\left(4\right)}$ and $E_{\left(6\right)}$ up to the eleventh digit is identical to $E_{\left(5\right)}$. The approximation $E_{\left(6\right)}$, also, up to the seventh digit is identical to $E_{\left(3\right)}$, given by (118). However, if the approximate value for $E_{\left(3\right)}$ is compared with the value for $E_{\left(6\right)}$, the coincidence only exists up to the fifth digit.

It is, therefore, necessary to check what the impact is of the approximations in the eccentric anomaly E on the approximations of the time of propagation of the signal.

8.2. Numerical Calculation of the First $O\left({\displaystyle \frac{1}{c}}\right)$ Correctionin the Propagation Time

The corresponding elliptic integral (59) of the second kind for the eccentric anomaly approximation $E_{\left(3\right)}$ is

(122)$T_1^{\left(E_3\right)}=\underset{0}{\overset{E_{\left(3\right)}}{\int }}\sqrt[.]{1-e^2{cos}^{2}E}dE=-0.317557268125933936045.$

If we compare this value with the value (118) $E_{\left(3\right)}=-0.31758547588467897473$ [rad], then it can be noted that the integration changes the value of $E_{\left(3\right)}$ after the fourth digit after the decimal dot.

The same integral with $E_{\left(6\right)}$ as an upper integration limit is

(123)$T_1^{\left(E_6\right)}=\underset{0}{\overset{E_{\left(6\right)}}{\int }}\sqrt[.]{1-e^2{cos}^{2}E}dE=-0.317558568963886638536.$

This value is different from the value (122) after the fifth digit, so it is a better approximation. However, if compared with $E_{\left(6\right)}=$$-0.317585484114454995929$ [rad], the integration in (123) changes the value of $E_{\left(6\right)}$ after the fourth digit after the decimal dot.

The calculation of the first $O\left({\displaystyle \frac{1}{c}}\right)$ time correction ${\displaystyle \frac{a}{c}}{T}_1$ for the values of the $GPS$ orbit and for the eccentric anomaly $E_{\left(3\right)}$ gives the following numerical value:

(124)${\displaystyle \frac{a}{c}}{T}_1^{\left(E_3\right)}=-0.0281341332790419\left[\mathrm{s}\right].$

Correspondingly, the first $O\left({\displaystyle \frac{1}{c}}\right)$ time correction ${\displaystyle \frac{a}{c}}{T}_1^{\left(E_6\right)}$ for the eccentric anomaly $E_{\left(6\right)}$ is

(125)${\displaystyle \frac{a}{c}}{T}_1^{\left(E_6\right)}=-0.0281342485273829\left[\mathrm{s}\right].$

So, the two $O\left({\displaystyle \frac{1}{c}}\right)$ time corrections (124) and (125) are identical up to the sixth digit after the decimal dot. Consequently, the important conclusion is that the eccentric anomaly at the sixth approximation level can ensure microsecond stability of the $O\left({\displaystyle \frac{1}{c}}\right)$ time correction. This means also that the third approximation $E_{\left(3\right)}$ can be reliably used, at least at the microsecond level.

One more argument in support of the microsecond approximation is that for a satellite, moving along circular orbit (when $e=0$), the above expression acquires the following form:

${\displaystyle \frac{a}{c}}{T}_{1\left(circular\right)}^{\left(E_3\right)}={\displaystyle \frac{a}{c}}{E}_{\left(3\right)}=-0.028136517830159266\left[\mathrm{s}\right]$

(126)$=-28136.517830159266\left[\mathsf{\mu }\mathrm{s}\right].$

Evidently, the coordinate time for the circular orbit differs from the coordinate time for the $GPS$ elliptic orbit only after the fifth digit after the decimal dot. This is, however, a clear proof that in order to account for changes in the propagation time in the micro- ($\mathsf{\mu }$) and the nano-range, the ellipticity of the orbit should be accounted in the calculation of the propagation time.

8.3. Numerical Calculation of the Shapiro Delay Term (the $O\left({\displaystyle \frac{1}{c^3}}\right)$ Time Correction)

Now let us calculate numerically the whole under-integral expression in formulae (41) for the propagation time, taking into account the $25-$th digit after the decimal dot

(127)${\displaystyle \frac{a}{c}}\sqrt{1-e^2{cos}^{2}E}-{\displaystyle \frac{2GM}{c^3}}\sqrt[.]{{\displaystyle \frac{1+ecos{E}_{\left(3\right)}}{1-ecos{E}_{\left(3\right)}}}}$

(128)$=-0.0885884561709019072629880\left[\mathrm{s}\right]-0.0000000000299618094618168\left[\mathrm{s}\right]$

(129)$=-0.0885884562008637167248048\left[\mathrm{s}\right].$

The second Shapiro delay term contains ten zeroes after the decimal dot, so the overall result of the calculation of both terms in (127) is changed by the Shapiro term at and after the nanosecond level (1 ns = ${10}^{-9}$ s). The second Shapiro delay term in (128) is equal to $2.99618\times {10}^{-9}\left[\mathrm{s}\right]\approx 2.99618\left[\mathrm{ns}\right]$, which is smaller than the timing error for a receiver at the geometric center of a tetrahedron satellite configuration, measured to be 10 [ns] [49]. Here it should be clarified that the calculation in (128) and in (129) is of an approximate value because one should take the numerical values for T after performing the integration of the elliptic integral. For example, from formulae (125)

${\displaystyle \frac{a}{c}}{T}_1^{\left(E_6\right)}=-0.0281342485273829\left[\mathrm{s}\right],$

obtained after integrating the first term in (127), it is seen that the first number $0.0885884561709019072629880$[s] in (128) has been diminished to the number $0.0281342485273829\left[\mathrm{s}\right]$ in (125), which represents a substantial change. In the literature (see [50]) time transfer has been discussed between a geostationary satellite and a station on the rotating Earth. This time transfer amounts to $80\left[\mathrm{ps}\right]=0.008\left[\mathrm{ns}\right]$ for the $c^{-3}$ terms, accounting for the Shapiro time delay (in [50] called a gravitational time delay). However, two-way time transfer between two stations at Braunschweig (Germany) and FTZ (Darmstadt, Germany) can be not so accurate and can amount to $300\left[\mathrm{ps}\right]=0,3\left[\mathrm{ns}\right]$, the uncertainties due to the indeterminacies in the position of the emitter (the station) and the receiver (the satellite), or vice versa. From the numerical calculation for the above case of inter-satellite communications, the relativistic time transfer (i.e., the propagation time of the signal) will be of another magnitude, depending also on the given final value for $E^{\left(6\right)}$. For the moment, we do not have a numerical estimation of the second integral term in (41), consequently we take the whole under-integral expression.

It should be noted also that in the previous calculations of the iterations $E_{\left(3\right)}$ and $E_{\left(6\right)}$ of the eccentric anomaly angle, it was taken with a negative sign, due to the accepted convention. After obtaining the result of the calculation of $T_1^{\left(E_3\right)}$ and $T_1^{\left(E_6\right)}$ from the integrals (122) and (123), they were with a negative sign. This means that after the integrating of the first term in (127), there will be a negative sign before the first number in (128) and the negative sign of the second term in (129) remains. This is so because after calculating ${\displaystyle \frac{2GM}{c^3}}\sqrt[.]{{\displaystyle \frac{1+ecos{E}_{\left(3\right)}}{1-ecos{E}_{\left(3\right)}}}}$, the obtained number is positive. Thus, the two numbers in (128) should be added, giving a total number of $-0.0885884562008637167248048$ in (129), greater than the first number in (128) after the ninth digit after the decimal dot.

Although approximate, this result, based on the simple application of the celestial mechanics approach in the calculation of the propagation time, is interesting, if compared with the calculations, based on the formulae (1) $T={\displaystyle \frac{R_{AB}}{c}}+2{\displaystyle \frac{G_{\oplus }{M}_E}{c^3}}ln\left({\displaystyle \frac{r_A+r_B+R_{AB}}{r_A+r_B-R_{AB}}}\right)$. From the literature, it can be seen that relativistic effects on light propagation from satellite laser ranging (SLR) data are measured with an accuracy of 1 $\mathsf{\mu }\mathrm{m}$ (1 micrometer), which corresponds to $0.01$ $\mathrm{p}\mathrm{s}$ (see also paper [50]). The relativistic observables, defined for the $GRAIL$ mission in [6], are also accurate to 1 $\mathsf{\mu }\mathrm{m}$. No doubt that the numerical calculation of the integral over the Shapiro delay term in (128) will give a more realistic result.

8.4. Propagation Time of the Signal Compared to the Celestial Time—Consistency of the Numerical Calculation

Let us compare the calculated propagation time of the signal with the celestial time $t_{cel}$, which can be calculated from $M=n(t-t_P)$ (M is the mean anomaly, given in the numerical data), where $t_P$ is the initial time of perigee passage. The idea is to calculate the celestial time of motion for the given value M of the mean anomaly and, thus, to compare the distance travelled by the satellite with the distance travelled by the light (laser or radio) signal. Since this calculation, also, will be approximate, we have taken only the value $T_1^{\left(E_6\right)}$ for the propagation time from the integral (123) and (125):

$T_1^{\left(E_6\right)}={\displaystyle \frac{a}{c}}\underset{0}{\overset{E_{\left(6\right)}}{\int }}\sqrt[.]{1-e^2{cos}^{2}E}dE=-0.0281342485273829\left[\mathrm{s}\right],$

which does not account for the second term, responsible for the Shapiro delay. Let us take the ratio of the propagation time to the celestial time:

(130)${\displaystyle \frac{\mathrm{propagation}\mathrm{time}\mathrm{for}\mathrm{electromagnetic}\mathrm{signals}}{\mathrm{celestial}\mathrm{time}\mathrm{from}\mathrm{the}\mathrm{Kepler}\mathrm{equation}}}={\displaystyle \frac{0.0281342485273829\left[\mathrm{s}\right]}{37.508256148\left[\mathrm{s}\right]}}$

(131)$=0.0007500814864964887.$

In other words, the propagation time for the signal is about ${10}^4$ times smaller than the celestial time, which will mean that for the celestial time 37.5082561 [s], the satellite moves 145,125 [km] (the velocity of the satellite is taken to be $v_S=3.874$$[\mathrm{km}/\mathrm{s}]$) and the light signal will move at a distance of

(132)$26,558.1510169176263505735093268096\left[\mathrm{km}\right].$

This value is obtained after multiplying the light velocity $c=299,792.458$$\left[{\displaystyle \frac{\mathrm{km}}{\mathrm{s}}}\right]$ with the calculated propagation time, taking into account the sixth iteration $T_1^{\left(E_6\right)}$ according to (125).

It is interesting to compare this value for the distance travelled by the signal with the value, obtained after taking only the first term in (128), for the propagation time, where the term ${\displaystyle \frac{a}{c}}\sqrt{1-e^2{cos}^{2}E}$ has not been integrated. This term does not account for the second term for the Shapiro delay, which contained 10 zeroes after the decimal dot and was of the order $0.299$ ns (nanoseconds). The result will be

(133)$26,558.151025899950855259224944504\left[\mathrm{km}\right].$

This number is greater (after the fourth digit after the decimal dot) than the preceding number (132), so the integration of the expression (125)

${\displaystyle \frac{a}{c}}\underset{0}{\overset{E_{\left(6\right)}}{\int }}\sqrt[.]{1-e^2{cos}^{2}E}dE$

in fact diminishes the propagation time (only the first $O\left({\displaystyle \frac{1}{c}}\right)$ correction) in comparison with the first term ${\displaystyle \frac{a}{c}}\sqrt{1-e^2{cos}^{2}E}$ in the expression (128). The difference between the numbers (133) and (132) is

(134)$0.0000089823245046857156176944\left[\mathrm{km}\right]\approx 0.898\left[\mathrm{cm}\right]\approx 8.98\left[\mathrm{mm}\right].$

In both cases of the numerical calculation of the numbers (132) and (133), the numerical result confirms that the distance travelled by the light signal (numbers (132) and (133)) is much greater than the distance of $145.125\left[\mathrm{km}\right]$ travelled by the satellite. Respectively, for the calculated propagation time $0.0281342485273829$$\left[\mathrm{s}\right]$ (only the first $O\left({\displaystyle \frac{1}{c}}\right)$ correction $T_1^{\left(6\right)}$ (125)), the satellite will move at a distance of

(135)${\displaystyle \frac{145.125\left[\mathrm{km}\right]}{37.5082561\left[\mathrm{s}\right]}}\times 0.0281342485273829\left[\mathrm{s}\right]$

(136)$=0.1088555758671073849925\left[\mathrm{km}\right].$

This distance is considerably smaller than the distance travelled by the signal. In fact, it is about 108.8 m, which is $0.000004098763343794755826290507\sim 0.409\times {10}^{-5}$ times smaller than the distance, given by the number (132)

$26,558.151016917626350573509326809\left[\mathrm{km}\right],$

travelled by the light or radio signal for the calculated propagation time of 0.0281342485273829 s.

8.5. Comparison of the Result in Formulae (134) witha Result in Other Papers. Approximations in the Calculations of the ShapiroDelay Term

The result of 8.98 [mm] in (134) was obtained by substracting the propagation distance, obtained from the first term in the formula (41), initially, by taking into account the integrated term, and afterwards, by calculating the numerical value only for the under-integral expression (see also (128)).

However, in the known review article by Ashby [51,52], it was proved that at the level of several millimeters, spatial curvature effects have to be considered, which is evident from the simple equality

(137)$\int dr\left[1+{\displaystyle \frac{G{M}_E}{c^{2}r}}\right]\approx r_2-r_1+{\displaystyle \frac{G{M}_E}{c^2}}ln{\displaystyle \frac{r_2}{r_1}}.$

In this paper, the second term was estimated to be $4.43\times ln\left(4.2\right)\approx 6.3$ mm. It might be concluded, yet, that in the first term in (134), spatial curvature effects are important and have to be included. In the Shapiro formulae, the first term is the Euclidean distance, divided by the light velocity. From a formal point of view, it might be believed that the first terms ${\displaystyle \frac{a}{c}}\sqrt{1-e^2{cos}^{2}E}$ in (128) and ${\displaystyle \frac{a}{c}}\underset{0}{\overset{E_{\left(6\right)}}{\int }}\sqrt[.]{1-e^2{cos}^{2}E}dE$ in (125) should correspond to the first term ${\displaystyle \frac{R}{c}}$ in the Shapiro formulae. Therefore, since it might be suspected that curvature effects may appear yet in these formulaes, it is not quite correct to consider any correspondence ${\displaystyle \frac{a}{c}}\sqrt{1-e^2{cos}^{2}E}\Rightarrow {\displaystyle \frac{R}{c}}$ or ${\displaystyle \frac{a}{c}}\underset{0}{\overset{E_{\left(6\right)}}{\int }}\sqrt[.]{1-e^2{cos}^{2}E}dE\Rightarrow {\displaystyle \frac{R}{c}}$. In view of this, in a recent paper [25], a new modification of the Shapiro formulae has been proposed, where the first term ${\displaystyle \frac{R_{.}}{c}}$ is replaced by the term ${\displaystyle \frac{1}{c}}\underset{path}{\int }\sqrt{{\left(\stackrel{.}{x}\right)}^2+{\left(\stackrel{.}{y}\right)}^2+{\left(\stackrel{.}{z}\right)}^2}ds$. In this term, s is some parameter, which parametrizes some curve, for example, the space curve, determining the space coordinates on the satellite orbit.

8.6. Comparison of the Propagation Time with the Atomic Time of theAtomic Clocks of the Satellites in the Near-Earth Space—Consistency of theResults

Now, let us point out one another numerical fact, which is interesting, because it confirms the consistency of the approach of calculating the propagation time on the base of using the parametrization of the satellite orbit. In the book [53] it was pointed out that at an altitude of 20,184 km, due to the difference in the gravitational potential, the atomic time of the satellite clock runs faster by 45 $\mathsf{\mu }\mathrm{s}/\mathrm{d}$ (microseconds per day, 1 $\mathsf{\mu }\mathrm{s}={10}^{-6}$ s). Consequently, for one second, the atomic time will run faster by $0.5208333\times {10}^{-9}$ s. In the preceding sections, it was proved that the dominant part in the propagation time $dT$ (i.e., the part without the Shapiro delay term, which is very small) is $dT=0.0281341332790419$ $\left[\mathrm{s}\right]$, given by the number (125) in Section 8.2. Corresponding to this propagation time interval $dT$, the atomic time is $d\tau =0.0146531934\times {10}^{-9}$ s and can easily be found from 45 $\mathsf{\mu }\mathrm{s}/\mathrm{d}$ (keeping in mind how many seconds are contained in 24 h, i.e., 1 day). The ratio of the two time intervals is equal to

(138)${\displaystyle \frac{d\tau }{dT}}={\displaystyle \frac{0.0146531934.{10}^{-9}}{0.028134248527382}}=0.520831163687866092.{10}^{-9}.$

The very small atomic time interval, compared to the propagation interval, means that the atomic time can serve as a standard for measuring the propagation time, because it will be able to detect changes even at the nanosecond level.

It should be noted that this value is of the order of ${10}^{-9}$, but it is a little greater than the value $\beta ={\displaystyle \frac{2V}{c^2}}=0.334\times {10}^{-9}$.

Similarly, if a clock on an orbit is compared to a clock on the Earth’s surface, due to the net effect of time dilation and gravitational redshift, the satellite clock will appear to run fast by 38 $\mathsf{\mu }\mathrm{s}/\mathrm{d}$ [53], which also is considered to be an enormous difference for an atomic clock with precision of a few nanoseconds. So, for this value, a similar number will be obtained for the ratio ${\displaystyle \frac{d\tau }{dT}}$, again confirming that this nanosecond $({10}^{-9}$$\mathrm{s})$ precision of the atomic clock (no matter whether it is on the Earths surface or on a satellite in a near-Earth orbit) is sufficient for measuring the propagation time T of the signal.

9. Justification of the Theoretical Approach in ThisPaper from the Viewpoint of General Relativity Theory and DifferentialGeometry

As mentioned in the introduction, the basic approach of the calculation of the propagation time in this paper is based on finding the solution for the propagation time from the null cone equation and parametrizing the space coordinates for both cases, considered as plane elliptical orbits or space-oriented orbits by the orbital coordinates. So, in fact, the light-like geodesics are not used in this approach. In the previous chapter, it was demonstrated that the solution for the propagation time does not give any time, related to the motion of the satellite, but the calculated propagation time is of the order of microseconds (1 $\mathsf{\mu }\mathrm{s}={10}^{-6}$ s) or nanoseconds ($1\mathrm{ns}={10}^{-9}$ s). This is the first fact in favor of the correctness of using the approach of the null cone equation. The second fact is related to the proper dimension of seconds of the expression for the propagation time for the case of a signal, emitted by a satellite on a plane elliptical orbit.

Now we shall present the third argument, related to a theorem, proved in the book on general relativity theory by Fock [19]. This theorem is related to known ones from basic textbooks on differential geometry [24,26]. Let us first begin with the theorem in [19]

Theorem 1. 

Let $F={\displaystyle \frac{1}{2}}{g}_{\alpha \beta }{\displaystyle \frac{d{x}^{\alpha }}{ds}}$${\displaystyle \frac{d{x}^{\beta }}{ds}}$. $F=0$ is the null cone equation and the extremum of the integral $s=\underset{s_1}{\overset{s_2}{\int }}Lds$ is searched $(L=\sqrt{2F})$; then, the Lagrange equation of motion will give

(139) ${\displaystyle \frac{dF}{ds}}=0⟹F=const⟹{\displaystyle \frac{d}{ds}}{\displaystyle \frac{\partial F}{\partial {\stackrel{.}{x}}_{\alpha }}}-{\displaystyle \frac{\partial F}{\partial x_{\alpha }}}=0.$

From the last equation, the geodesic equations $(\alpha ,\beta ,\nu =0,1,2,3)$

(140) ${\displaystyle \frac{d^2{x}_{\nu }}{d{s}^2}}+{\Gamma }_{\alpha \beta }^{\nu }{\displaystyle \frac{d{x}_{\alpha }}{ds}}{\displaystyle \frac{d{x}_{\beta }}{ds}}=0$

can be obtained. Since the constant in (139) can be set up to zero, this would mean that the geodesic Equation (140) is compatible with the null cone equation $F={\displaystyle \frac{1}{2}}{g}_{\alpha \beta }{\displaystyle \frac{d{x}^{\alpha }}{ds}}$${\displaystyle \frac{d{x}^{\beta }}{ds}}=0$.

It is important to acquire the proper understanding of this statement: It means that every solution of the null equation (including in the proper parametrization of the differentials, related to the trajectory of motion of the satellite) is also a solution of the light-like geodesic Equation (140). In the “backward direction”, however, the statement is not true, meaning that there might be other solutions of the geodesic equation, which are not solutions of the null cone equation.

The meaning of the geodesic as the shortest path, joining two given points (but only in the local sense, i.e a sufficiently short part of the geodesics), is confirmed also by the following theorem in the monograph by Prasolov [26]. Let us present the exact formulation of this theorem:

Theorem 2. 

Any point p on a surface S has a neighborhood U such that the length of a geodesic γ going from the point p to some point q and lying entirely in the neighborhood U does not exceed the length of any curve α on S, joining the points p and q. Moreover, if the lengths of the curves γ and α are equal, then γ and α coincide as non-parametrizable curves.

In fact, what does it mean that “the length of the geodesics from point p to point q does not exceed the length of any other curve” (but only locally, i.e., in the vicinity of the given neighborhood)? This means that this geodesic line is minimal. So, this confirms the theorem in another, third well-known book on differential geometry by Mishtenko and Fomenko [24]:

Theorem 3. 

The geodesic line $\gamma :\left[a,b\right]\to M^n$ is minimal if it is not longer than any smooth path, joining the endpoints $\gamma \left(a\right)$ and $\gamma \left(b\right)$.

So, if we have, for example, the functional $E\left(\gamma \left(s\right)\right)=\int \mid {\stackrel{.}{\gamma }}^2\left(s\right)\mid ds=\int g_{ij}\left(x\right){\stackrel{.}{x}}^i\left(s\right){\stackrel{.}{x}}^j\left(s\right)ds$ for the null cone equation, and we are searching for the extremals, these extremals will be derived from the geodesics—the multitude of the solutions for the extremals from the geodesics turns out to be larger in comparison with the one for the other local curves in the vicinity, in particular, for the functional $E\left(\gamma \right(s\left)\right)$ for the null cone equation.

On the basis of combining the results in the previously mentioned monographs by Fock [19,26], one easily establishes the validity of the following theorem in [24]:

Theorem 4. 

The extremals of the length functional $L\left(\gamma \left(s\right)\right)=\int \sqrt[.]{g_{ij}\left(x\right){\stackrel{.}{x}}^i\left(s\right){\stackrel{.}{x}}^j\left(s\right)}ds$ are smooth functions, derived from the geodesics (by means of smooth change of the parameter along the curve). Then, each extremal of the functional $E\left(\gamma \left(s\right)\right)=\int \mid {\stackrel{.}{\gamma }}^2\left(s\right)\mid ds=\int g_{ij}\left(x\right){\stackrel{.}{x}}^i\left(s\right){\stackrel{.}{x}}^j\left(s\right)ds$ (called also action functional of the trajectory) is also an extremal of the functional $L\left(\gamma \right(s\left)\right)=$$\int \sqrt[.]{g_{ij}\left(x\right){\stackrel{.}{x}}^i\left(s\right){\stackrel{.}{x}}^j\left(s\right)}ds$. However, each extremal of the functional $L\left(\gamma \right(s\left)\right)$ is not necessarily an extremal of the functional $E\left(\gamma \right(s\left)\right)$.

There is a simple inequality $L^2\left(\gamma \left(s\right)\right)\le$ $E\left(\gamma \right(s\left)\right)$ [24] in confirmation of the above theorem, which can be immediately proved by means of the Schwartz inequality

(141)$\left(\underset{0}{\overset{1}{\int }}f\left(s\right)g\left(s\right)ds\right)\le \left(\underset{0}{\overset{1}{\int }}{f}^2\left(s\right)ds\right)\left(\underset{0}{\overset{1}{\int }}{g}^2\left(s\right)ds\right)$

for the following functions [24]:

(142)$f\left(s\right)\equiv 1g\left(s\right)\equiv \mid \stackrel{.}{\gamma }\left(s\right)\mid .$

It should be noted that the parameter s in the book [19] is defined simply as a “space parameter along a curved line”. In the above formulaes, s is also a space parameter. But as noted in [24], when $g\left(s\right)=const$, i.e., $\gamma \left(s\right)=const$ and the parameter s is proportional to the length of the arc of the curved line, an equality can be derived $L^2\left(\gamma \left(s\right)\right)=E\left(\gamma \left(s\right)\right)$. Thus, every solution of the null cone equation $g_{ij}\left(x\right){\stackrel{.}{x}}^i\left(s\right){\stackrel{.}{x}}^j\left(s\right)=0$ (i.e., a zero-extremal of $E\left(\gamma \right(s\left)\right)$) will also be a zero extremal of $L^2\left(\gamma \left(s\right)\right)$. If this zero extremal corresponds to the propagation of the signal along the geodesics of the minimal length, the propagation time will satisfy Fermat’s principle regarding the least time of propagation of light between two space points. Consequently, even if there are other extremals for the length functional, of interest are only those, related to the minimal geodesics, obtained from the action functional $E\left(\gamma \left(s\right)\right)=\int g_{ij}\left(x\right){\stackrel{.}{x}}^i\left(s\right){\stackrel{.}{x}}^j\left(s\right)ds$, which ensures the fulfillment of Fermat’s principle.

10. New Approach for Analytical Calculation of Elliptic Integrals—Analytic Relation Between Elliptic Integrals in the Weierstrass and in the Legendre Form

In this section, some new analytical algorithms for the calculation of elliptic integrals will be presented, which have been published in the paper [33]. The main advantage of the theoretical method is that (1) it will be applicable for various values of the modulus q of the elliptic integral, no matter whether q is a small or a big number, and (2) the integrals, encountered in the calculation of the propagation time in the preceding sections, were in the Legendre form, which contains the specific fourth-order polynomial $(1-y^2)(1-q^2{y}^2)$ under the square root in the denominator.

However, in the previous Section 4, the definitions for elliptic integrals with arbitrary polynomials of the third and of the fourth degree were given. Such integrals are very often encountered in cosmology, nonlinear evolution equations, black hole physics, etc. For example, in the definition for the luminosity in cosmology, the more general polynomial of the fourth degree under the square root in the denominator was used. Also, there is an interesting interplay between the methods of the standard three-dimensional uniformization and of the four-dimensional uniformization, described in the monograph by Whittaker and Watson [35]. This method will be exposed in detail further and will allow not only the analytical computation of the Weierstrass invariants $g_2$ and $g_3$ for four-dimensional elliptic integrals of the type $J_n^{\left(4\right)}\left(y\right)$ but also the comparison of elliptic integrals in the Weierstrass form $J_n^{\left(3\right)}\left(x\right)$, both of them defined by formulaes (48)–(51) in Section 4.1, concerning the definitions of higher-order elliptic integrals.

Some physical applications of elliptic integrals were presented in the paper [33].

10.1. Transforming an Elliptic Integral in the Legendre Form into anElliptic Integral in the Weierstrass Form

The purpose of this section will be to propose a transformation, which will transform the elliptic integral

(143)${\tilde{J}}_0^{\left(4\right)}\left(y\right)=\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}$

in the Legendre form into an elliptic integral in the Weierstrass form

(144)${\tilde{J}}_0^{\left(3\right)}\left(x\right)=-\sqrt{a}\int {\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}.$

For the purpose, the transformation

(145)$x={\displaystyle \frac{a}{y^2}}+b$

shall be applied. It is possible to calculate the parameters a and b, so that the integral (144) is satisfied. These parameters represent complicated functions of the modulus parameter q of the integral

(146)$a={\displaystyle \frac{9g_3}{g_2}}K\left(q\right),$

(147)$b=-{\displaystyle \frac{a}{3}}(1+q^2)=-{\displaystyle \frac{3g_3}{g_2}}K\left(q\right).$

In (147), $K\left(q\right)$ is a rational function, depending on the modulus parameter q of the elliptic integral in the Legendre form

(148)$K\left(q\right)\equiv {\displaystyle \frac{q^4-q^2+1}{2q^4-5q^2+2}}.$

However, $g_2$ and $g_3$ remain undetermined. This will turn out to be an advantage of the applied formalism, since, further, another representation of the integral in the Weierstrass form shall be found with different Weierstrass invariants ${\overline{g}}_2$ and ${\overline{g}}_3$, which will depend explicitly on the modulus parameter q. Then, by comparing both representations in the Weierstrass form and proving a new theorem, $g_2$ and $g_3$ will be expressed through ${\overline{g}}_2$ and ${\overline{g}}_3$ in a complicated way.

Further, the following important property will be very important.

Corollary 1. 

The parameter b satisfies the cubic equation, i.e.,

(149) $4b^3-g_{2}b-g_3=0.$

The proof is straightforward, if the expressions for $g_2$ and $g_3$ are expressed from (146) and (147) and are substituted into the above cubic equation.

10.2. Equivalence Between a Formulae from Integral Calculus and the Recurrent System of Equations for the Elliptic Integral in Terms of the y-Variable

Due to the property (149), the elliptic integral (144) can be represented in the equivalent way:

(150)${\tilde{J}}_0^{\left(3\right)}\left(\overline{x}\right)=-{\displaystyle \frac{\sqrt{a}}{2}}\int {\displaystyle \frac{d\overline{x}}{{\overline{x}}^{\frac{1}{2}}\sqrt{{\overline{x}}^2+3b\overline{x}+(b^2-\frac{g_2}{4})}}},$

where $\overline{x}\equiv x-b$. Such integrals in the monograph by Timofeev [46] are analytically computed according to the formulae

$\int {\displaystyle \frac{d\overline{x}}{{\overline{x}}^m\sqrt{{\left(\tilde{a}+2\tilde{b}\overline{x}+\tilde{c}{\overline{x}}^2\right)}^{2n+1}}}}=-{\displaystyle \frac{\sqrt{{\left(\tilde{a}+2\tilde{b}\overline{x}+\tilde{c}{\overline{x}}^2\right)}^{2n-1}}}{(m-1)\tilde{a}{\overline{x}}^{m-1}}}$

$-{\displaystyle \frac{(2m+2n-3)\tilde{b}}{(m-1)\tilde{a}}}\int {\displaystyle \frac{d\overline{x}}{{\overline{x}}^{m-1}\sqrt{{\left(\tilde{a}+2\tilde{b}\overline{x}+\tilde{c}{\overline{x}}^2\right)}^{2n+1}}}}$

(151)$-{\displaystyle \frac{(m+2n-2)\tilde{c}}{(m-1)\tilde{a}}}\int {\displaystyle \frac{d\overline{x}}{{\overline{x}}^{m-2}\sqrt{{\left(\tilde{a}+2\tilde{b}\overline{x}+\tilde{c}{\overline{x}}^2\right)}^{2n+1}}}}.$

Since this is a formulae from integral calculus, it is interesting to prove the following theorem, which most surprisingly establishes the equivalence between the above relation from integral calculus with some of the recurrent system of equations for the case of a fourth-degree polynomial in the denominator of the elliptic integral [34]. The recurrent system for the case of a cubic polynomial under the square root has been investigated in the monographs of Fichtenholz [38] and Smirnov [45]. Below, the formulation of the newly proved theorem is given.

Theorem 5. 

The integral (151) is equivalent to the second equation

(152) ${\overline{z}}^{.}\tilde{Z}=3\tilde{c}{J}_4^{\left(4\right)}\left(\overline{z}\right)+4\tilde{b}{J}_2^{\left(4\right)}\left(\overline{z}\right)+\tilde{a}{J}_0^{\left(4\right)}\left(\overline{z}\right).$

from the recurrent system of equations

(153) ${\overline{z}}^m\tilde{Z}=(m+2)\tilde{c}{J}_{m+3}^{\left(4\right)}\left(\overline{z}\right)+(m+1)2\tilde{b}{J}_{m+1}^{\left(4\right)}\left(\overline{z}\right)+m\tilde{a}{J}_{m-1}^{\left(4\right)}\left(\overline{z}\right),$

where the fourth-degree polynomial $\tilde{Z}$ is defined as

(154) ${\tilde{Z}}^2=a_0{\overline{z}}^4+4a_1{\overline{z}}^3+6a_2{\overline{z}}^2+4a_3\overline{z}+a_4.$

The coefficients in the above polynomial have the values

(155) $a_0=\tilde{c},a_1=0,6a_2=2\tilde{b},a_3=0,a_4=\tilde{a}.$

The corresponding elliptic integrals of the first kind and of the third kind and of the n-th order are

(156) $J_n^{\left(4\right)}\left(\overline{z}\right)=\int {\displaystyle \frac{{\overline{z}}^n}{\sqrt{\tilde{a}+2\tilde{b}{\overline{z}}^2+\tilde{c}{\overline{z}}^{4.}}}}d\overline{z}=\int {\displaystyle \frac{{\overline{z}}^{n}d\overline{z}}{\sqrt[.]{{\overline{z}}^4+3b{\overline{z}}^2+(3b^2-\frac{g_2}{4})}}},$

(157) $H_n^{\left(4\right)}\left(\overline{z}\right)=\int {\displaystyle \frac{d\overline{z}}{{\overline{z}}^n\sqrt{\tilde{a}+2\tilde{b}{\overline{z}}^2+\tilde{c}{\overline{z}}^{4.}}}}=\int {\displaystyle \frac{d\overline{z}}{{\overline{z}}^n\sqrt[.]{{\overline{z}}^4+3b{\overline{z}}^2+(3b^2-\frac{g_2}{4})}}}.$

Due to the established equivalence, it can be concluded that the integral (151) cannot be used for the calculation of the zero-order elliptic integral in the Legendre form (143), because the recurrent relation (152) contains also the higher-order elliptic integrals $J_2^{\left(4\right)}\left(\overline{z}\right)$ and $J_4^{\left(4\right)}\left(\overline{z}\right)$.

10.3. A New Theorem, Concerning the Equivalence Between Some of theEquations in the Recurrent System

Theorem 6. 

From (153), the corresponding equation for $m=-3$

(158) ${\displaystyle \frac{\tilde{y}}{y^3}}=-q^2{J}_0\left(y\right)+2(1+q^2)H_2^{\left(4\right)}\left(y\right)-3H_4^{\left(4\right)}\left(y\right)$

can be transformed into the equation for $m=1$

(159) $y{\tilde{y}}^{\left(q\right)}=J_0\left(y\right)-2(1+q^2)J_2\left(y\right)+3q^2{J}_4\left(y\right)$

by means of the transformation $y\Rightarrow {\displaystyle \frac{1}{qy}}$ and the simple relations (for the concrete case for $m=2$)

(160) $H_m\left(y\right)=-q^m{J}_m\left({\displaystyle \frac{1}{qy}}\right),H_m\left({\displaystyle \frac{1}{qy}}\right)=-q^m{J}_m\left(y\right).$

In such a way, only one independent equation can be obtained from the $4D$ recurrent system of equations

(161) $2(1+q^2)J_2^{\left(q\right)}\left(y\right)=3q^2{J}_4^{\left(q\right)}\left(y\right)+J_0^{\left(q\right)}\left(y\right)-y{\tilde{y}}^{\left(q\right)}.$

10.4. Relations Between the Recurrent System of Equations in Terms of the y and x Variables

10.4.1. The Recurrent System of Equations in Terms of the x Variable

Now, we shall write the recurrent system of equations for the elliptic integrals

$I_m^{\left(3\right)}\left(x\right)\equiv \int {\displaystyle \frac{x^{m}dx}{\sqrt{4x^3-g_{2}x-g_3}}},$

(162)$G_m^{\left(3\right)}\left(x\right)\equiv \int {\displaystyle \frac{dx}{x^m\sqrt{4x^3-g_{2}x-g_3}}},x={\displaystyle \frac{{\tilde{a}}_2}{y^2}}+{\tilde{b}}_2,$

as this has been performed in the monographs [38,45]. Generally, the polynomial of a fourth-order under the square root is denoted as

(163)$Y^2=a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4,$

but for the concrete case of the integrals (162), the coefficients will be given by

(164)$a_0=0,a_1=1,a_2=0,a_3=-{\displaystyle \frac{g_2}{4}},a_4=-g_3.$

For this choice of the coefficients, due to $a_0=0$, the polynomial becomes a cubic one, so we change the variables in (163) from $(Y,y)$ to $(X,x)$, i.e., $X^2=4x^3-g_{2}x-g_3$, where the Weierstrass invariants (coefficients) $g_2$ and $g_3$ will be defined in the next sections.

The recurrent system of equations is easily found to be (omitting the upper indice “$\left(3\right)$” in $I_m^{\left(3\right)}$)

(165)$x^{n}X=2(2n+3)I_{n+2}-{\displaystyle \frac{g_2}{2}}(2n+1)I_n-n{g}_3{I}_{n-1}.$

It should be noted also that the variable transformation $x={\displaystyle \frac{{\tilde{a}}_2}{y^2}}+{\tilde{b}}_2$ in (162) contains other parameters ${\tilde{a}}_2$ and ${\tilde{b}}_2$, which are different from the parameters $a_2$ and $b_2$ in the transformation (145). Now, we shall transform the three-degree polynomials under the square root in (162) into the known four-degree polynomials for the elliptic integral in the Legendre form. The transformation $x={\displaystyle \frac{{\tilde{a}}_2}{y^2}}+{\tilde{b}}_2$ (as well as the previously applied transformation (145) $x={\displaystyle \frac{a}{y^2}}+b$) are generalizations of the transformation $x=e_3+{\displaystyle \frac{(e_1-e_3)}{y^2}}$ in the paper [54], where $e_1,e_2,e_3$ are constant roots of a cubic polynomial. The transformation transforms the integral $\int {\displaystyle \frac{dx}{\sqrt{4(x-e_1)(x-e_2)(x-e_3)}}}$ into the integral $-{\displaystyle \frac{1}{\sqrt{e_1-e_3}}}\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}$, where $q^2={\displaystyle \frac{(e_2-e_3)}{(e_1-e_3)}}$. Under another transformation, $y\Rightarrow {\displaystyle \frac{1}{\tilde{k}y}}$, the x-variable will transform as

(166)$\tilde{x}\equiv {\tilde{a}}_2{\tilde{k}}^2{y}^2+{\tilde{b}}_2,$

where a new variable $\tilde{x}$ has been introduced. The two recurrent systems of equations, (153) (written in terms of the y-variable) and the recurrent system (165) for the x-variable, are not independent because

(167)$I_0^{\left(3\right)}\left(x\right)=-{\displaystyle \frac{J_0^{\left(q\right)}\left(y\right)}{{\left(\sqrt{a}\right)}^q}}={\displaystyle \frac{2{\tilde{a}}_2}{{\left(\sqrt{a}\right)}^q}}{J}_0^{\left(\tilde{k}\right)}\left(y\right),$

where $\tilde{k}$ is determined so that $4x^3-g_2\left(q\right)x-g_3\left(q\right)$ is transformed into the expression $(1-y^2)(1-$${\tilde{k}}^2{y}^2)$. Consequently, the last expression is satisfied if $\tilde{k}$ is defined as

$\tilde{k}=-1-{\displaystyle \frac{6\alpha }{\sqrt[3]{2}}},\alpha ={\tilde{b}}_2\mathrm{is}\mathrm{a}\mathrm{root}\mathrm{of}4{\tilde{b}}_2^3-g_2^{\left(q\right)}{\tilde{b}}_2-g_3^{\left(q\right)}=0,$

(168)${\tilde{a}}_2={\displaystyle \frac{1}{\sqrt[3]{4}}}={\displaystyle \frac{{\tilde{k}}^4-{\tilde{k}}^2+1}{3g_2^{\left(q\right)}}}.$

Two central problems remain to be solved:

How to combine the derived equations in terms of the x variable with the found relation (161);

How to estimate the parameter $g_2$ in Equation (168). In fact, further, in the next sections, it will be proved that $g_2=g_2^{\left(q\right)}\equiv g_2\left(q\right)$ and $g_3=g_3^{\left(q\right)}\equiv g_3\left(q\right)$ will depend in a complicated way on the modulus parameter q, but this dependence will be found from another Weierstrass representation of the integral $J_0^{\left(4\right)}\left(y\right)$ (143). This new representation will not contain the conformal factor $\sqrt{a}$ in Equation (144), but instead will be characterized by different Weierstrass invariants ${\overline{g}}_2^{\left(q\right)}$ and ${\overline{g}}_3^{\left(q\right)}$, which will be explicitly determined by the method of the s.c. “fourth-degree polynomial uniformization”. This algebraic method will be considered in detail in the following sections.

10.4.2. Combining the Two Recurrent Systems of Equations

From the two recurrent system of equations in terms of the y-variable and the x-variable, the following equation can be obtained:

(169)${\tilde{X}}^{\left(\tilde{k}\right)}=N_1({\tilde{a}}_2,{\tilde{b}}_2,\tilde{k})J_2^{\left(\tilde{k}\right)}\left(y\right)+N_2({\tilde{a}}_2,{\tilde{b}}_2,\tilde{k})J_0^{\left(\tilde{k}\right)}\left(y\right)+N_3({\tilde{a}}_2,{\tilde{b}}_2,\tilde{k})=$

(170)$=N_1^{\left(\tilde{k}\right)}{J}_2^{\left(\tilde{k}\right)}\left(y\right)+N_2^{\left(\tilde{k}\right)}\left[-{\displaystyle \frac{1}{2{\tilde{a}}_2}}{J}_0^{\left(q\right)}\left(y\right)\right]+N_3^{\left(\tilde{k}\right)},$

where

(171)$N_1({\tilde{a}}_2,{\tilde{b}}_2,\tilde{k})\equiv N_1^{\left(\tilde{k}\right)}\equiv 8{\tilde{a}}_2^3{\tilde{k}}^2(1+{\tilde{k}}^2)+24{\tilde{a}}_2^2{\tilde{b}}_2{\tilde{k}}^2,$

(172)$N_2({\tilde{a}}_2,{\tilde{b}}_2,\tilde{k},)\equiv N_2^{\left(\tilde{k}\right)}\equiv 12{\tilde{a}}_2{\tilde{b}}_2^2-{\displaystyle \frac{g_2{\tilde{a}}_2}{{\left(\sqrt{a}\right)}^{\left(q\right)}}}-4{\tilde{a}}_2^3{\tilde{k}}^2,$

(173)$N_3({\tilde{a}}_2,{\tilde{b}}_2,\tilde{k})\equiv N_3^{\left(\tilde{k}\right)}\equiv 4{\tilde{a}}_2^3{\tilde{k}}^2{y}^2{\tilde{y}}^{\left(\tilde{k}\right)}.$

In (169), the expression for ${\tilde{X}}^{\left(\tilde{k}\right)}$ is defined as ${\tilde{X}}^{\left(\tilde{k}\right)}\equiv 4{\tilde{x}}^3-g_2^{\left(\tilde{k}\right)}\tilde{x}-g_3^{\left(\tilde{k}\right)}$. The notations $g_2^{\left(\tilde{k}\right)}$ and $g_3^{\left(\tilde{k}\right)}$ mean that these Weierstrass invariants are expressed as complicated functions of the parameter $\tilde{k}$.

Now, it can be noted that if it is possible to make the replacement $\tilde{k}⟺q$, then another equation will be obtained:

(174)${\tilde{X}}^{\left(q\right)}=N_1^{\left(q\right)}{J}_2^{\left(q\right)}\left(y\right)+N_2^{\left(q\right)}{J}_0^{\left(q\right)}\left(y\right)+N_3^{\left(q\right)},$

and in this way, the Equations (170) and (174) will give the following solution for the zero-order elliptic integral in the Legendre form:

$J_0^{\left(4\right)}\left(y\right)\equiv J_0^{\left(q\right)}\left(y\right)={\displaystyle \frac{2{\tilde{X}}^{\left(q\right)}}{\sqrt[3]{2}{N}_2^{\left(q\right)}}}-{\displaystyle \frac{N_3^{\left(q\right)}}{N_2^{\left(q\right)}}}$

(175)$-{\displaystyle \frac{\sqrt[3]{2}{\tilde{X}}^{\left(\tilde{k}\right)}\left(\sqrt[3]{2}-2N_2^{\left(q\right)}\right)}{N_2^{\left(q\right)}\left(\sqrt[3]{2}-2\right)}}-{\displaystyle \frac{N_1^{\left(q\right)}\left(2{\tilde{X}}^{\left(q\right)}-\sqrt[3]{2}{N}_3^{\left(q\right)}\right)}{N_2^{\left(q\right)}\sqrt[3]{2}}}.$

As a consequence of the replacement $\tilde{k}⟺q$, one also has the following equality, which follows from the defining relations (168):

(176)${\displaystyle \frac{{\tilde{k}}^4-{\tilde{k}}^2+1}{3g_2\left(\tilde{k}\right)}}={\displaystyle \frac{q^4-q^2+1}{3g_2\left(q\right)}}.$

If this is a quite a general condition (not imposing any restrictions), then the replacement $\tilde{k}⟺q$ will be justified. Evidently, this will depend on the determination of $g_2^{\left(q\right)}\equiv g_2\left(q\right)$, which will be given in the next sections.

10.5. Another Integral in the Weierstrass Form After Applying theFour-Dimensional Uniformization

Further, by “four-dimensional uniformization”, it shall be meant that the polynomial of the fourth degree under the square root in the denominator of the elliptic integral shall be “uniformized” by introducing complicated functions, depending on the Weierstrass function and its derivative, which represent functions of the complex variable z. The term “uniformization” shall be clarified in the next section.

10.6. The Standard Method for $3D$ Uniformization in EllipticFunctions Theory—Reminder of the Basic Facts

In accord with the standard notions in elliptic functions theory [34,55], the s.c. “fundamental parallelogram” is defined on the two-periodic lattice (see also the monograph [56])

(177)$\Lambda \equiv \left\{m{\omega }_1+n{\omega }_2;0\le m\le 1,0\le n\le 1,Im\left({\displaystyle \frac{{\omega }_1}{{\omega }_2}}\right)>0\right\}.$

The two periods ${\omega }_1$ and ${\omega }_2$ constitute the basis on the complex plane and, thus, each complex number can be represented as $z=$ $m{\omega }_1+n{\omega }_2$. Then, the Weierstrass function

(178)$\rho \left(z\right)\equiv {\displaystyle \frac{1}{z^2}}+\sum \left({\displaystyle \frac{1}{{(z-\omega )}^2}}-{\displaystyle \frac{1}{{\omega }^2}}\right)$

is an example of an elliptic function. By definition, $f\left(z\right)$ is an elliptic function if it is meromorphic and doubly periodic. The double periodicity of the elliptic function means that the equality

(179)$f(z+n{\omega }_1+m{\omega }_2)=f\left(z\right)$

is fulfilled, where n and m are integer numbers. A complex function is meromorphic, if it does not have singular points other than poles [34]. The singular points can be of three kinds [57,58]:

Removable singular points, for which ${lim}_{z\to a}f\left(z\right)$ exists and is a finite one. These points constitute the s.c. main part in the Laurent decomposition

(180)$f\left(z\right)=c_0+c_1(z-a)+.......+c_n{(z-a)}^n\Rightarrow \underset{z\to a}{lim}f\left(z\right)=c_0.$

Essentially singular points, for which ${lim}_{z\to a}f\left(z\right)$ does not exist.

The case with poles, for which ${lim}_{z\to a}f\left(z\right)=\infty$. For such a case, the Laurent decomposition is

(181)$f\left(z\right)={\displaystyle \frac{c_{-n}}{{(z-a)}^n}}+........+{\displaystyle \frac{c_{-1}}{z-a}}+\sum _{k=0}^{\infty }{c}_k{(z-a)}^k.$

The first derivative of the Weierstrass function $\rho \left(z\right)$ (178) is

(182)${\rho }^{\prime }\left(z\right)\equiv -{\displaystyle \frac{2}{z^3}}-2\sum {\displaystyle \frac{1}{{(z-\omega )}^3}}$

and it satisfies the parametrizable form of the cubic algebraic equation

(183)$y^2=4x^3-g_{2}x-g_3\mathrm{i}.\mathrm{e}.,x\equiv \rho \left(z\right),y\equiv {\rho }^{\prime }\left(z\right),$

where the functions $\rho \left(z\right)$ and ${\rho }^{\prime }\left(z\right)$ are called also uniformization functions for the above cubic algebraic curve $y^2=4x^3-g_{2}x-g_3$.

The summation in the defining formulaes (178) for the Weierstrass function $\rho \left(z\right)$ is over all points on the two-dimensional lattice on the complex plane, i.e., on all possible numbers m and n. The coefficient functions $g_2$ and $g_3$ are the s.c. “invariants of the Weierstrass function” $\rho \left(z\right)$

(184)$g_2=60G_4=60\sum {\displaystyle \frac{{}_{ }{}^{\prime }1_{ }^{ }}{{\omega }^4}},g_3=140G_6=140\sum {\displaystyle \frac{{}_{ }{}^{\prime }1_{ }^{ }}{{\omega }^6}},$

defined as infinite convergent sums over the period $\omega$ of the two-periodic lattice. The prime “′” above the two sums means that the period $\omega =0$ is excluded from the summation, so that the expressions would not tend to infinity.

10.7. Application of the Weierstrass Integral and of the WeierstrassElliptic Curve in the Parametrizable Form

This section has the purpose of clarifying the importance of finding the Weierstrass invariants $g_2$ and $g_3$ for the calculation of an elliptic integral in the Weierstrass form $\int {\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}$. In the next sections on the basis of the methods for four-dimensional uniformization, a primary goal will be to find the Weierstrass invariants $g_2$ and $g_3$ not only for integrals in the Weierstrass form, but also for integrals in the Legendre form and also for more general integrals (50) (for $n=0$)

$J_0^{\left(4\right)}\left(y\right)=\int {\displaystyle \frac{dy}{\sqrt{a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4}}},$

with an arbitrary polynomial of the fourth-degree under the square in the denominator. Here, it shall be clarified that finding the Weierstrass invariants will be possible after transforming the integral initially given in the Legendre form into an integral in the Weierstrass form, and in such a way, the Weierstrass invariants $g_2$ and $g_3$ will depend in a complicated manner on the modulus q of the integral in the Legendre form. Thus, as it will be demonstrated by a theorem in the monograph by [59], for known $g_2$ and $g_3$, it will be possible to determine the solution (and the roots) of the cubic equation $y^2=4x^3-g_{2}x-g_3$. Then, the periods ${\omega }_1$ and ${\omega }_2$ of the two-dimensional lattice of the fundamental parallelogram (177) can be calculated.

Let us first note that finding a solution $z-z_0$ of the integral $\int {\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}$ in the Weierstrass form is equivalent to knowing the roots of the cubic polynomial $y^2=4x^3-g_{2}x-g_3$. This “equivalency” can be written as [57]

(185)$z-z_0=\underset{x_0}{\overset{x}{\int }}{\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}\iff y^2=4x^3-g_{2}x-g_3.$

In fact, the equality $x=\rho \left(z\right)$ is a solution also of the elliptic integral in (185) and this is the known problem of “inversion” for elliptic integrals. In (183), $\rho \left(z\right)$ is the Weierstrass function [57]:

(186)$\rho \left(z\right)\equiv {\displaystyle \frac{1}{z^2}}+\sum \left({\displaystyle \frac{1}{{(z-\omega )}^2}}-{\displaystyle \frac{1}{{\omega }^2}}\right)={\displaystyle \frac{1}{z^2}}+{\displaystyle \frac{g_2}{20}}{z}^2+{\displaystyle \frac{g_3}{28}}{z}^4+...$

From the second representation (186) for the Weierstrass function $\rho \left(z\right)$, it becomes clear why, in the literature, it is denoted by $\rho (z;g_2,g_3)$. The following theorem, proved in Ch. 6 of the monograph by Knapp [59], allows us to understand how the result of integration for the elliptic integral (143) ${\tilde{J}}_0^{\left(q\right)}\left(y\right)=\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}$ can be represented in an analytical form by means of the Weierstrass integral, after calculating the Weierstrass invariants $g_2$ and $g_3$ by means of the method of four-dimensional uniformization. This method will be presented in the next sections.

10.8. A Theorem about the Unique Correspondence Between the Periods ${\omega }_1$ and ${\omega }_2$ on the Two-DimensionalComplex Lattice and the Weierstrass Invariants $(g_2,g_3)$ for EllipticIntegrals in the Weierstrass Form

Theorem 7 

([59]). There exists an unique correspondence $(g_2(\Lambda ),g_3(\Lambda ))\iff \Lambda$ between the lattices Λ on the complex plane $\mathbb{C}$ and the pair $(g_2,g_3)$ of the complex numbers, the Weierstrass invariants, such that the discriminant $\Delta \equiv g_2^3-27g_3^2$ of the cubic polynomial $4x^3-g_{2}x-g_3=0$ is different from zero. If $e_1,e_2,e_3$ are the roots of the above polynomial, i.e.,

(187)$4x^3-g_{2}x-g_3=4(z-e_1)(z-e_2)(z-e_3),$

then the periods $({\omega }_1,{\omega }_2)$ on the complex lattice Λ (177) can be found from the following integrals:

(188)${\omega }_1=\underset{{\Gamma }_1}{\int }{\displaystyle \frac{dz}{2\sqrt{(z-e_1)(z-e_2)(z-e_3)}}},$

(189)${\omega }_2=\underset{\Gamma 2}{\int }{\displaystyle \frac{dz}{2\sqrt{(z-e_1)(z-e_2)(z-e_3)}}}.$

In the first integral, the unique branch ${\Gamma }_1$ of the square root is chosen with cuts from $e_1$ to $e_2$ and from $e_3$ to ∞, and in the second integral, the branch ${\Gamma }_2$ is with cuts from $e_2$ to $e_3$ and from $e_1$ to ∞.

Moreover, for every lattice $\Lambda$ on the complex plane, defined according to (177), a biholomorphic mapping is determined as $\phi :\mathbb{C}/\Lambda \mapsto E\left(\mathbb{C}\right)$, where $E\left(\mathbb{C}\right)$ denotes the elliptic curve $y^2=4x^3-g_2(\Lambda )x-g_3(\Lambda )$. The mapping $\phi$ is defined by means of the Weierstrass elliptic function $\rho \left(z\right)$ and its derivative ${\rho }^{\prime }\left(z\right)$ and the inverse mapping—by the corresponding elliptic integral.

Since $\rho \left(z\right)$ and ${\rho }^{\prime }\left(z\right)$ are uniformization functions, they will satisfy the cubic polynomial [34]

(190)${\left({\rho }^{\prime }\left(z\right)\right)}^2=4\left(\rho \left(z\right)-e_1\right)\left(\rho \left(z\right)-e_2\right)\left(\rho \left(z\right)-e_3\right),$

from where the following relations are easily derived [34], following from the algebraic properties of the roots of a cubic equation

(191)$e_1+e_2+e_3=0,e_1{e}_2+e_2{e}_3+e_1{e}_3=-{\displaystyle \frac{g_2}{4}},e_1{e}_2{e}_3={\displaystyle \frac{g_3}{4}}.$

Now, it may easily be calculated that

(192)$g_2^3-27g_3^2=16{(e_1-e_2)}^2{(e_2-e_3)}^2{(e_1-e_3)}^2.$

If the discriminant $\Delta \equiv g_2^3-27g_3^2$ (a well known notion from higher algebra) is different from zero, then the cubic equation $4x^3-g_2(\Lambda )x-g_3(\Lambda )=0$ will not have coinciding roots, i.e.,

(193)$e_1\ne e_2\ne e_3.$

This property is very important also for the theory of elliptic curves (see also the monograph [60]) since if the property $g_2^3-27g_3^2\ne 0$ is fulfilled, then a lattice $\Lambda$ on the complex plane exists, so that the Weierstrass invariants $g_2$ and $g_3$ can be defined as (184)

$g_2=60G_4=60\sum {\displaystyle \frac{{}_{ }{}^{\prime }1_{ }^{ }}{{\omega }^4}},g_3=140G_6=140\sum {\displaystyle \frac{{}_{ }{}^{\prime }1_{ }^{ }}{{\omega }^6}}.$

At this point, it becomes clear why the correspondence $(g_2(\Lambda ),g_3(\Lambda ))\iff \Lambda$ in the cited monograph [59] is in both directions—previously, we proved that if the lattice $\Lambda$ is given, then the Weierstrass invariants $g_2$ and $g_3$ can be defined, i.e., $\Lambda \Rightarrow (g_2(\Lambda ),g_3(\Lambda ))$.

10.8.1. The Consistency Problem for a Non-Zero Discriminant for the Case, When the Invariants in the Weierstrass Integrals are Expressed Through the Modulus Parameter q of the Integral in the Legendre Form

Now, the statement is proved in the opposite directions—if the Weierstrass invariants are given, then a complex lattice can be defined $(g_2(\Lambda ),g_3(\Lambda ))\Rightarrow \Lambda$. Note that this statement is not trivial and it can be reformulated in the following way: If $g_2$ and $g_3$ are given (real or complex) numbers or even functions of some parameter (such examples will be shown further), then a two-dimensional lattice $\Lambda$ on the complex plane with periods ${\omega }_1$ and ${\omega }_2$ can be defined, so that the sums in the defining equalities for $g_2$ and $g_3$ in (184) will be convergent. This is also a problem, which needs further investigation.

Further, it shall be understood why this is important in view of the fact that an elliptic integral in the Legendre form (with a modulus parameter q) shall be transformed to an elliptic integral in the Weierstrass form $\int {\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}$. The Weierstrass integral is characterized by invariants $g_2$ and $g_3$, which shall be expressed by complicated polynomials, depending on the modulus parameter q, i.e., $g_2=g_2\left(q\right)$ and $g_3=g_3\left(q\right)$. So, a natural question arises: Will the integral in the Weierstrass form with the calculated values for the Weierstrass invariants $g_2$ and $g_3$ be correctly defined? The answer will be affirmative, because from the above correspondence, it will follow that the lattice $\Lambda$ can be defined. However, if $g_2$ and $g_3$ are defined as $g_2=g_2\left(q\right)$ and $g_3=g_3\left(q\right)$, then the condition $\Delta \equiv g_2^3-27g_3^2\ne 0$ for the discriminant should also be fulfilled, which, in fact, will mean that the higher-order $g_2^3-27g_3^2$ polynomial of q should have no zeroes! Further, in the section for four-dimensional uniformization, we shall obtain the parametrization for $g_2$ and $g_3$ (212) in the form of the following higher-order polynomials of the modulus parameter q of the elliptic integral (143) in the Legendre form

${\overline{g}}_2^{(1,2)}\equiv q^2+{\displaystyle \frac{1}{12}}{(1+q^2)}^2>0$

and (213) (the notation in both equalities will be ${\overline{g}}_2$ and ${\overline{g}}_3$)

${\overline{g}}_3^{(1,2)}\equiv {\displaystyle \frac{q^6}{4}}-{\displaystyle \frac{q^4}{6}}-{\displaystyle \frac{5q^2}{12}}-{\displaystyle \frac{1}{6^3}}{\left[5q^2-1\right]}^3.$

Then, will the higher-order polynomial ${\overline{g}}_2^{3(1,2)}-27{\overline{g}}_3^{2(1,2)}$ have the property to be without any roots? The same question appears with respect to the found second pair of Weierstrass invariants (226) (the same as (212)),

${\overline{g}}_2=q^2+{\displaystyle \frac{1}{12}}{(1+q^2)}^2>0$

and also (227)

${\overline{g}}_3=-{\displaystyle \frac{1}{6}}{q}^2(1+q^2)+{\displaystyle \frac{{(1+q^2)}^3}{6^3}}.$

In view of the theorems, which have been proved, it follows that in both cases, the higher-order polynomial ${\overline{g}}_2^{3(1,2)}-27{\overline{g}}_3^{2(1,2)}$ should have no roots! Since both the polynomials ${\overline{g}}_2^{3(1,2)}-27{\overline{g}}_3^{2(1,2)}$ depend on higher degrees of q, which is defined for $0\le q<1$, it is known from higher algebra that there are many theorems [61] investigating the case when an arbitrary polynomial of the n-th order will have (or will not have) roots in the circle $q<1$. One such theorem, formulated in the monograph by Obreshkoff [61], is the theorem by Schur [62], published in 1918. The concrete application of this theorem requires cumbersome analytic calculations, because it should be applied to a chain of algebraic equations of diminishing degrees. For example, for the case of the polynomials (226) and (227), the polynomial ${\overline{g}}_2^{3(1,2)}-27{\overline{g}}_3^{2(1,2)}$ will be of the 12th degree! For polynomials of the fourth degree, the Schur theorem has been applied several times in the papers [20,40] and the results have shown excellent consistency with the physical meaning, related to these polynomials. Moreover, in [20], it has been demonstrated that the Schur theorem can be used both for the case of proving that a certain polynomial has roots in the circle $q<1$ and for the case of proving that it does not have such roots. In the next sections, we shall give only the formulation of this theorem.

On the other hand, if in one of the cases for the polynomial ${\overline{g}}_2^{3(1,2)}-27{\overline{g}}_3^{2(1,2)}$, it is proved that there are roots, then this will mean that either the corresponding parametrization should be rejected, or that the elliptic integral in the Legendre form is not defined for the values of q, satisfying the above higher-order polynomial ${\overline{g}}_2^{3(1,2)}-27{\overline{g}}_3^{2(1,2)}=0$. In standard monographs on elliptic functions and integrals, this interesting problem is not investigated.

Now, it remains to prove that the three roots $e_1,e_2,e_3$ are different and they will correspond to three different values for the Weierstrass function $\rho \left(z_1\right)$, $\rho \left(z_2\right)$, $\rho (z_1+z_2)$, corresponding to the values

(194)$z_1={\displaystyle \frac{{\omega }_1}{2}},z_2={\displaystyle \frac{{\omega }_2}{2}},z_1+z_2={\displaystyle \frac{{\omega }_1}{2}}+{\displaystyle \frac{{\omega }_2}{2}}.$

The proof is standard and can be found in the known monographs by [34,57,60].

A detailed study on whether the periods $({\omega }_1,{\omega }_2)$ are real and imaginary depending on the roots of the polynomial $4x^3-g_{2}x-g_3=0$ and on the properties of the integrals (188) and (189) is given in Ch. $6.5$ of the known monograph [63].

10.8.2. A Differential Equation for the Weierstrass Function as aConsequence of the Group-Theoretical Law for Summing up Points on the CubicCurve

Since the group law for summing up points on the cubic curve and the resulting differential equation will be applied further for finding other formulaes for the Weierstrass invariants $g_2$ and $g_3$ in (228), we shall review briefly the derivation of this equation, given in [34].

Let $e_1,e_2,e_3$ denote the three roots of the cubic curve $y^2=4x^3-g_{2}x-g_3$. If a straight line $y=fx+h$ ($f,h$ are coefficients) intersects the cubic curve, then, from the cubic equation

(195)${(fx+h)}^2=4x^3-g_{2}x-g_3,$

also from the Wiets formulae for the roots and from the equalities (195) for the three values of the Weierstrass function at the points $z=z_1$, $z=z_2$ and $z=z_3$, it follows [34] that

(196)$\rho \left(z_1\right)+\rho \left(z_2\right)+\rho (z_1+z_2)={\displaystyle \frac{f^2}{4}}.$

Since $x=\rho \left(z\right)$ and $y={\rho }^{\prime }\left(z\right)$ are uniformization functions for the cubic curve and also for the straight line $y=fx+h$, let us write it in terms of the Weierstrass function

(197)${\rho }^{\prime }\left(z_1\right)=f\rho \left(z_1\right)+h,{\rho }^{\prime }\left(z_2\right)=f\rho \left(z_2\right)+h.$

Expressing, from here, the coefficient f as

(198)$f={\displaystyle \frac{{\rho }^{\prime }\left(z_1\right)-{\rho }^{\prime }\left(z_2\right)}{\rho \left(z_1\right)-\rho \left(z_2\right)}}$

and substituting into the cubic Equation (195), the following differential equation is obtained, which is a consequence from the group-theoretic law for summing up points on a cubic curve:

(199)$\rho (z_1+z_2)=\rho \left(z_1\right)+\rho \left(z_2\right)+{\displaystyle \frac{1}{4}}{\left({\displaystyle \frac{{\rho }^{\prime }\left(z_1\right)-{\rho }^{\prime }\left(z_2\right)}{\rho \left(z_1\right)-\rho \left(z_2\right)}}\right)}^2.$

Closely related with the group theoretical law on elliptic curves is the problem of rational points on elliptic curves. Both topics are extensively discussed in the monographs [34,64,65].

11. Uniformization of Fourth-Degree Algebraic Equations andApplication to the Theory of Elliptic Integrals—New Representation of anIntegral in the Legendre Form as an Integral in the Weierstrass Form

The derivation of the second representation of the elliptic integral in the Weierstrass form is based on the following theorem, proved in the monograph by Whittaker and Watson [35].

Theorem 8. 

Let

(200) $J_0^{\left(4\right)}\left(y\right)=\int {\displaystyle \frac{dy}{\sqrt{a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4}}}$

be a zero-order elliptic integral with the fourth-degree polynomial under the square root in the denominator

(201) $Y\equiv f\left(y\right)\equiv a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4,f\left(y_0\right)=0.$

Then, after the variable changes

(202) ${\displaystyle \frac{1}{t-y_0}}={\displaystyle \frac{\sigma -\frac{1}{2}{A}_2}{A_3}},{\displaystyle \frac{1}{y-y_0}}={\displaystyle \frac{s-\frac{1}{2}{A}_2}{A_3}}$

the integral can be rewritten as

(203) $I_0^{\left(3\right)}\left(x\right)={\int }_s^{\infty }{\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-{\overline{g}}_2\sigma -{\overline{g}}_3}}},$

where the Weierstrass invariants ${\overline{g}}_2$ and ${\overline{g}}_3$ can be exactly calculated as

(204) ${\overline{g}}_2=3A_2^2-4A_1{A}_3;{\overline{g}}_3=2A_1{A}_2{A}_3-A_2^3-A_0{A}_3^2$

and $A_0,A_1,A_2,A_3$ are the introduced notations for the polynomial expressions

(205) $A_0\equiv a_0,A_1\equiv a_0{y}_0+a_1,$

(206) $A_2\equiv a_0{y}_0^2+2a_1{y}_0+a_2,A_3\equiv a_0{y}_0^3+3a_1{y}_0^2+3a_2{y}_0+a_3.$

It can easily be verified that if the expressions (205) and (206) for $A_0,A_1,A_2,A_3$ are used, then the invariant ${\overline{g}}_2$ in (204) can be calculated as

(207)${\overline{g}}_2=a_0{a}_4-4a_1{a}_3+3a_2^2.$

Later, in the other formulation of the theorem in the next paragraph, this formulae will be used, but it will be taken from the theorem in the monograph [35]. Note that under the changes in the coefficients

(208)$4a_1\to a_1,6a_2\to a_2,4a_3\to a_3,a_4\to a_4$

in the fourth-degree polynomial (201), the invariant ${\overline{g}}_2$ changes as

(209)${\overline{g}}_2=a_0{a}_4-{\displaystyle \frac{a_1{a}_3}{4}}+{\displaystyle \frac{a_2^2}{12}}.$

Technically, the calculation of the invariant ${\overline{g}}_3$ is more difficult.

11.1. Four-Dimensional Uniformization of the (Zero-Order) EllipticIntegral in the Legendre Form and its Representation as an Integral in theWeierstrass Form

First, we shall clarify that the application of the above formulaes (202)–(206) to the elliptic integral (200) with a polynomial of the fourth degree in the square root and transforming it to the elliptic integral $I_0^{\left(3\right)}\left(x\right)$ (203) will be called a “four-dimensional uniformization”. If the fourth-degree polynomial (201)

$Y\equiv f\left(y\right)\equiv a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4,f\left(y_0\right)=0.$

is represented in a parametrizable form by functions y and Y, expressed in a complicated way by functions, depending on the Weierstrass function $\rho \left(z\right)$ and its derivative ${\rho }^{\prime }\left(z\right)$, then this will be called “an uniformization of the fourth-degree polynomial (201) with the uniformization functions $y=y(y_0,z,\rho \left(z\right),{\rho }^{\prime }\left(z\right))$ and $Y=Y(y_0,z,\rho \left(z\right),{\rho }^{\prime }\left(z\right))$”. Further, these uniformization functions will be shown.

In the concrete case for the polynomial in the Legendre form

(210)$f\left(y\right)\equiv (1-y^2)(1-q^2{y}^2)=1-(1+q^2)y^2+y^4,$

one can compute ${\overline{g}}_2^{(1,2)}$ and ${\overline{g}}_3^{(1,2)}$ for the roots of the polynomial $f\left(y\right)$ (201) $y_{1,2}=\pm 1$ and $y_{3,4}=\pm {\displaystyle \frac{1}{q}}$ and also for the coefficients $a_0,a_1,a_2,a_3,a_4$, taking into account the identity of the polynomials (201) and (210):

(211)$a_0\equiv a_1\equiv a_3\equiv 0,a_2\equiv -{\displaystyle \frac{1}{6}}(1+q^2),a_4\equiv 1.$

The two pairs of coefficient functions ${\overline{g}}_2^{(1,2)}$ and ${\overline{g}}_3^{(1,2)}$, calculated according to the formulaes (204)–(206) and depending on functions of the modulus parameter q, are

(212)${\overline{g}}_2^{(1,2)}\equiv q^2+{\displaystyle \frac{1}{12}}{(1+q^2)}^2>0,$

(213)${\overline{g}}_3^{(1,2)}\equiv {\displaystyle \frac{q^6}{4}}-{\displaystyle \frac{q^4}{6}}-{\displaystyle \frac{5q^2}{12}}-{\displaystyle \frac{1}{6^3}}{\left[5q^2-1\right]}^3$

$=\left({\displaystyle \frac{1}{4}}-{\displaystyle \frac{5^3}{6^3}}\right)q^6+\left(-{\displaystyle \frac{1}{6}}+{\displaystyle \frac{75}{6^3}}\right)q^4$

(214)$-\left({\displaystyle \frac{1}{12}}+{\displaystyle \frac{15}{6^3}}\right)q^2+{\displaystyle \frac{1}{6^3}}.$

Similar expressions for ${\overline{g}}_2^{(3,4)}$ and ${\overline{g}}_3^{(3,4)}$ can also be written, but they will not be given here. Here the important fact, related to the theory of elliptic functions, is that if ${\overline{g}}_2^{(1,2)}$, ${\overline{g}}_3^{(1,2)}$ or ${\overline{g}}_2^{(3,4)}$, ${\overline{g}}_3^{(3,4)}$ are known, then, from the theorem in the monograph [59], it will follow that the Weierstrass function $\rho (z;{\overline{g}}_2,{\overline{g}}_3)$ is uniquely determined.

11.2. Another Formulation of the Theorem for the Uniformization ofAlgebraic Equations of the Fourth-Order and Application to the Theory ofElliptic Integrals

Theorem 9. 

If the Weierstrass elliptic function is defined as $\rho (z;{\overline{g}}_2,{\overline{g}}_3)$, then the two variables

(215) $y=y_0+{\displaystyle \frac{1}{4}}{\displaystyle \frac{f^{’}\left(y_0\right)}{\left[\rho (z;{\overline{g}}_2,{\overline{g}}_3)-\frac{1}{24}{f}^{’}\left(y_0\right)\right]}},$

(216) $Y=-{\displaystyle \frac{1}{4}}{\displaystyle \frac{f^{’}\left(y_0\right){\rho }^{’}(z;{\overline{g}}_2,{\overline{g}}_3)}{{\left[\rho (z;{\overline{g}}_2,{\overline{g}}_3)-\frac{1}{24}{f}^{’}\left(y_0\right)\right]}^2}}$

define an unique parametrization of the $4D$ algebraic curve (201)

$Y\equiv f\left(y\right)\equiv a_0{y}^4+4a_1{y}^3+6a_2{y}^2+4a_{3}y+a_4,f\left(y_0\right)=0.$

Consequently, the method for “four-dimensional uniformization” can be considered to be the analogue of the known “three-dimensional uniformization”, given by (183).

In the above formulaes, the point $y_0$ is supposed to be a root of the quartic polynomial (201), but this requirement can be relaxed and $y_0=a$ can be an arbitrary constant. The polynomial $f\left(y\right)$ is supposed not to have multiple roots. Then, the uniformization function y can be represented in the form given in the monograph [35] (the formulaes are shown as a part of two problems, after the formulation of the theorem for the two uniformization functions (215) and (216):

(217)$y=y_0+{\displaystyle \frac{F_1(y_0,z)}{F_2(y_0,z)}},$

where $F_1(y_0,z)$ and $F_2(y_0,z)$ are the expressions

$F_1(y_0,z)\equiv {\left[f\left(y_0\right)\right]}^{{\displaystyle \frac{1}{2}}}{\rho }^{’}(z;{\overline{g}}_2,{\overline{g}}_3)+{\displaystyle \frac{1}{2}}{f}^{’}\left(y_0\right)\{\rho (z;{\overline{g}}_2,{\overline{g}}_3)$

(218)$-{\displaystyle \frac{1}{24}}{f}^{’}\left(y_0\right)\}+{\displaystyle \frac{1}{24}}f\left(y_0\right){\displaystyle \frac{d^{3}f\left(y_{.}\right)}{d{y}^3}}{\mid }_{y=y_0},$

(219)$F_2(y_0,z)\equiv 2{\left[\rho (z;{\overline{g}}_2,{\overline{g}}_3)-{\displaystyle \frac{1}{24}}{f}^{’}\left(y_0\right)\right]}^2-{\displaystyle \frac{1}{48}}f\left(y_0\right)f^{IV}\left(y_0\right).$

If, again, $y_0$ is supposed to be a root of the algebraic equation $f\left(y\right)=0$ and the Weierstrass function $\rho (z;{\overline{g}}_2,{\overline{g}}_3)$ is known, then every other value of the polynomial $f\left(y\right)$ can be found by means of the following concise formulae ([35], also formulated as a separate problem):

(220)${\left[f\left(y\right)\right]}^{{\displaystyle \frac{1}{2}}}=-{\displaystyle \frac{f^{’}\left(y_0\right){\rho }^{’}(z;{\overline{g}}_2,{\overline{g}}_3)}{4\left[\rho (z;{\overline{g}}_2,{\overline{g}}_3)-\frac{1}{24}{f}^{’}\left(y_0\right)\right]}}.$

The formulaes (218)–(220) will not be used further in this paper, but probably, the interested reader can find some application of these expressions.

11.3. Second Formulation of the Theorem for Four-DimensionalUniformization

This theorem is also exposed in the monograph by E.T. Whittaker and G. N. Watson [35]. In another monograph by the same author Whittaker [66], many useful examples are given about the solutions in terms of elliptic functions of known equations from theoretical mechanics and analytical dynamics. We shall review the corresponding theorem by making a slight modification to the fourth-degree polynomial (201)—the numbers $4,6$ and 4 in the coefficients are removed.

Theorem 10. 

If $f\left(s\right)$ is a fourth-order polynomial

(221)$f\left(y\right)=a_0{y}^4+a_1{y}^3+a_2{y}^2+a_{3}y+a_4$

and the variable transformation $z=\underset{y_0}{\overset{y}{\int }}{\displaystyle \frac{ds}{\sqrt{f\left(s\right)}}}$ is performed in the $4D$ elliptic integral

(222)$J_0^{\left(4\right)}\left(y\right)=\int {\displaystyle \frac{dy}{\sqrt{a_0{y}^4+a_1{y}^3+a_2{y}^2+a_{3}y+a_4}}},$

then the Weierstrass invariants  ${\overline{g}}_2$  and  ${\overline{g}}_3$ in the integral (203)  $I_0^{\left(3\right)}\left(s\right)={\int }_s^{\infty }{\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-{\overline{g}}_2\sigma -{\overline{g}}_3}}}$  are uniquely calculated as

(223)${\overline{g}}_2=a_0{a}_4-{\displaystyle \frac{1}{4}}{a}_1{a}_3+{\displaystyle \frac{1}{12}}{a}_2^2,$

(224)${\overline{g}}_3={\displaystyle \frac{1}{6}}{a}_0{a}_2{a}_4+{\displaystyle \frac{1}{48}}{a}_1{a}_2{a}_3-{\displaystyle \frac{a_2^3}{6^3}}-{\displaystyle \frac{1}{16}}(a_0{a}_3+a_1^2{a}_4).$

One may note that in this formulation of the theorem, the Weierstrass invariants ${\overline{g}}_2$ and ${\overline{g}}_3$ do not depend on any root $y_0$ of the polynomial (221).

11.4. Four-Dimensional Uniformization of the Elliptic Integral in theLegendre Form and Its Representation in the Weierstrass Form—AnotherCalculation by Using the Second Formulation of the Theorem

Let us take, as an example, which will be used further, the polynomial

(225)$f\left(y\right)\equiv (1-y^2)(1-q^2{y}^2)=1-(1+q^2)y^2+y^4.$

The two Weierstrass invariants are calculated to be

(226)${\overline{g}}_2=q^2+{\displaystyle \frac{1}{12}}{(1+q^2)}^2>0,$

(227)${\overline{g}}_3=-{\displaystyle \frac{1}{6}}{q}^2(1+q^2)+{\displaystyle \frac{{(1+q^2)}^3}{6^3}}.$

Clearly, the result for ${\overline{g}}_2$ in inequality (226) coincides with the result for ${\overline{g}}_2^{(1,2)}\equiv q^2+{\displaystyle \frac{1}{12}}{(1+q^2)}^2>0$ in (212), but formulae (227) for ${\overline{g}}_3$ is quite different from formulaes (213) and (214) for ${\overline{g}}_3^{(1,2)}$. The same refers also to the formulae for ${\overline{g}}_3^{(3,4)}$, not presented here.

The important consequences are as follows:

For this case, there is only one pair of Weierstrass invariants ${\overline{g}}_2$ and ${\overline{g}}_3$ instead of the two pairs ${\overline{g}}_2^{(1,2)},{\overline{g}}_3^{(1,2)}$ and ${\overline{g}}_2^{(3,4)},{\overline{g}}_3^{(3,4)}$ for the previous case, considered in the preceding Section 11.

The expression for ${\overline{g}}_2$ (226) coincides with the expression for ${\overline{g}}_2^{(1,2)}$ in (212), but expression ${\overline{g}}_3$ in (227) is different from expression (213) for ${\overline{g}}_2^{(1,2)}$. It is not clear whether this will be so for the case of an elliptic integral with an arbitrary polynomial in the denominator.

The two invariants ${\overline{g}}_2$ and ${\overline{g}}_3$ are symmetric with respect to the interchanges $a_0\rightleftarrows a_4$ and also $a_1\rightleftarrows a_3$. The polynomial $f\left(y\right)$ is also invariant with respect to these changes.

12. Comparison of Elliptic Integrals in Different Variables and withDifferent Weierstrass Invariants

In this section, the method for calculating the Weierstrass invariants $g_2$ and $g_3$ will be presented, which is based on the combination of two approaches. In the first approach (see Section 10.1), the initial integral ${\tilde{J}}_0^{\left(q\right)}\left(y\right)$ was brought to the Weierstrass form (144) (multiplied by the conformal factor $\sqrt{a}=\sqrt{{\displaystyle \frac{9g_3}{g_2}}.{\displaystyle \frac{\left(q^4-q^2+1\right)}{\left(2q^4-5q^2+2\right)}}}$) by means of the transformation (145) $x={\displaystyle \frac{a}{y^2}}+b$, as a result of which, the integral acquires the form (144). This will be Equation (228) below.

The second approach is based on the method of “four-dimensional uniformization”, as a result of which, the integral is brought to the second representation (203) $I_0^{\left(3\right)}\left(x\right)={\int }_s^{\infty }{\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-{\overline{g}}_2\sigma -{\overline{g}}_3}}}$. This will be Equation (229) below.

As a result of the two representations, the following two equalities can be written:

(228)${\tilde{J}}_0^{\left(q\right)}\left(y\right)=\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}=-\sqrt{a}\int {\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}$

(229)$=\int {\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-{\overline{g}}_2\sigma -{\overline{g}}_3}}}.$

Now, we shall use the following theorem, proved in the monograph by Hancock [67], in order to transform the integral (228).

Theorem 11. 

The elliptic integral in the Weierstrass representation (144) (equation (A) in (228) in terms of the “x”-variable) can be written in terms of the variable σ as

(230)$\int {\displaystyle \frac{dx}{\sqrt{4x^3-g_{2}x-g_3}}}=-\int {\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-g_2\sigma -g_3}}}$

after performing the variable transformation

(231)$x=x_0+{\displaystyle \frac{(x_1-x_0)(x_2-x_0)}{\sigma -x_0}},$

where $x_0,x_1$ and $x_2$ are the roots of the cubic polynomial $4x^3-g_{2}x-g_3=0$.

In the original formulation of the theorem in [67] (page 29), the transformation $t=e_1+{\displaystyle \frac{(e_2-e_1)(e_3-e_1)}{s-e_1}}$ has been applied to the integral $\int {\displaystyle \frac{dt}{\sqrt{(t-e_1)}(t-e_2)(t-e_3)}}$. Note, however, the minus sign in front of the right term in Equation (230), which is important and can easily be confirmed by a direct calculation. Applying, again, the above theorem with respect to Equations (229) and also (203), the equality (230) can be written as

(232)$\int {\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-g_2\sigma -g_3}}}=-{\displaystyle \frac{1}{\sqrt[3]{a}}}\int {\displaystyle \frac{d\overline{\sigma }}{\sqrt{4{\overline{\sigma }}^3-{\widehat{g}}_2\overline{\sigma }-{\widehat{g}}_3}}},$

where the variable transformation $\sigma ={\displaystyle \frac{1}{\sqrt[3]{a}}}\overline{\sigma }$ has been performed and the notations

(233)${\widehat{g}}_2\equiv {\overline{g}}_2{a}^{{\displaystyle \frac{2}{3}}},{\widehat{g}}_3\equiv {\overline{g}}_{3}a$

are introduced. From Equations (230) and (232) (combined with (203)), we obtain

(234)$\int {\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-g_2\sigma -g_3}}}={\displaystyle \frac{1}{\sqrt[3]{a}}}\int {\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-{\widehat{g}}_2\sigma -{\widehat{g}}_3}}}.$

Note the important peculiarity of Equation (234): The Weierstrass invariants on the left-hand side of (234) are $g_2$, $g_3$, while on the right-hand side, they are ${\widehat{g}}_2\equiv {\overline{g}}_2{a}^{{\displaystyle \frac{2}{3}}},{\widehat{g}}_3\equiv {\overline{g}}_{3}a$.

12.1. Algebraic Roots of the Cubic Equations and Comparison of theWeierstrass Invariants—The Case of Positive Discriminant

From the Cardano formulae, we have, for the roots of the cubic polynomial $4{\sigma }^3-g_2\sigma -g_3=0$ ($i=1,2,3$), the following formulae (see the monographs [61,68]):

${\alpha }_1=A_1+B_1=\sqrt[3]{{\displaystyle \frac{g_3}{8}}+{\displaystyle \frac{1}{4}}\sqrt{{\displaystyle \frac{g_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{g}_2^3}}$

(235)$+\sqrt[3]{{\displaystyle \frac{g_3}{8}}-{\displaystyle \frac{1}{4}}\sqrt{{\displaystyle \frac{g_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{g}_2^3}}.$

This is only the first root of the cubic equation $4{\sigma }^3-g_2\sigma -g_3=0$. Further, we shall investigate the case of a positive discriminant $\Delta ={\displaystyle \frac{q^2}{4}}+{\displaystyle \frac{p^3}{27}}>0$ for the cubic equation $x^3+px+q=0$ [61]. For the given cubic equation $4{\sigma }^3-g_2\sigma -g_3=0$ with the first root (235), the discriminant is

(236)$\Delta ={\displaystyle \frac{g_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{g}_2^3={\displaystyle \frac{g_3^2}{4}}-{\displaystyle \frac{1}{108}}{g}_2^3>0\Rightarrow g_3^2>{\displaystyle \frac{1}{27}}{g}_2^3.$

Let us assume that the roots of this polynomial are ${\alpha }_1,{\beta }_1,{\gamma }_1$ and the roots of the polynomial $4{\sigma }^3-{\widehat{g}}_2\sigma -\widehat{g}=0$ are ${\alpha }_2,{\beta }_2,{\gamma }_2$. Then, the other two roots ${\beta }_1,{\gamma }_1$ can be found if the formulae for the s.c. n-th root of the number 1 is taken into account. This simple formulae is

(237)$\sqrt[n]{1}=\sqrt[n]{ϵ}=cos{\displaystyle \frac{2\pi }{n}}+isin{\displaystyle \frac{2\pi }{n}}.$

For $n=3$, and by means of the formulae

(238)$\sqrt[3]{1}=\sqrt[3]{ϵ}=cos{\displaystyle \frac{2\pi }{3}}+isin{\displaystyle \frac{2\pi }{3}}=-cos60+isin60=-{\displaystyle \frac{1}{2}}+{\displaystyle \frac{\sqrt[.]{3}}{2}},$

the roots ${\beta }_1,{\gamma }_1$ can be expressed as

${\beta }_1=A_1ϵ+B_1{ϵ}^2=A_1(-{\displaystyle \frac{1}{2}}+i{\displaystyle \frac{\sqrt[.]{3}}{2}})$

(239)$+B_1(-{\displaystyle \frac{1}{2}}-i{\displaystyle \frac{\sqrt[.]{3}}{2}})=-{\displaystyle \frac{(A_1+B_1)}{2}}+i{\displaystyle \frac{\sqrt[.]{3}}{2}}(A_1-B_1),$

${\gamma }_1=A_1{ϵ}^2+B_1{ϵ}^{.}=A_1(-{\displaystyle \frac{1}{2}}-i{\displaystyle \frac{\sqrt[.]{3}}{2}})$

(240)$+B_1(-{\displaystyle \frac{1}{2}}+i{\displaystyle \frac{\sqrt[.]{3}}{2}})=-{\displaystyle \frac{(A_1+B_1)}{2}}-i{\displaystyle \frac{\sqrt[.]{3}}{2}}(A_1-B_1).$

From the denominators of the two integrals in (234), we should have the equality of two cubic polynomials under a square root, the first expression for the polynomial obtained from the square of the second cubic polynomial, multiplied by $\sqrt[3]{a}$

(241)$\sqrt[.]{4{\sigma }^3-g_2\sigma -g_3}=\sqrt[3]{a}\sqrt[.]{4{\sigma }^3-{\widehat{g}}_2\sigma -{\widehat{g}}_3}.$

This expression can be written in the form

$(\sigma -{\alpha }_1)(\sigma -{\beta }_1)(\sigma -{\gamma }_1)={\left(\sqrt[3]{a}\right)}^2(\sigma -{\alpha }_2)(\sigma -{\beta }_2)(\sigma -{\gamma }_2)$

(242)$=(a^{{\displaystyle \frac{2}{9}}}\sigma -a^{{\displaystyle \frac{2}{9}}}{\alpha }_2)(a^{{\displaystyle \frac{2}{9}}}\sigma -a^{{\displaystyle \frac{2}{9}}}{\beta }_2)(a^{{\displaystyle \frac{2}{9}}}\sigma -a^{{\displaystyle \frac{2}{9}}}{\gamma }_2).$

The first root $a^{{\displaystyle \frac{2}{9}}}{\alpha }_2$ of the cubic equation on the right-hand side of (242), analogously to the Cardano formulae (235) for the previous equation, can be expressed as

$a^{{\displaystyle \frac{2}{9}}}{\alpha }_2=A_2+B_2=\sqrt[3]{{\displaystyle \frac{{\widehat{g}}_3}{8}}+{\displaystyle \frac{1}{4}}\sqrt{{\displaystyle \frac{{\widehat{g}}_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{\widehat{g}}_2^3}}$

(243)$+\sqrt[3]{{\displaystyle \frac{{\widehat{g}}_3}{8}}-{\displaystyle \frac{1}{4}}\sqrt{{\displaystyle \frac{{\widehat{g}}_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{\widehat{g}}_2^3}}.$

From here, after multiplication of the right-hand side of (243) with $a^{-{\displaystyle \frac{2}{9}}}$, only the root ${\alpha }_2$ can be expressed as

${\alpha }_2=a^{-{\displaystyle \frac{2}{9}}}(A_2+B_2)=\sqrt[3]{a^{-{\displaystyle \frac{2}{27}}}{\displaystyle \frac{{\widehat{g}}_3}{8}}+{\displaystyle \frac{1}{4}}\sqrt{a^{-{\displaystyle \frac{1}{27}}}{\displaystyle \frac{{\widehat{g}}_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{a}^{-{\displaystyle \frac{1}{27}}}{\widehat{g}}_2^3}}$

(244)$+\sqrt[3]{a^{-{\displaystyle \frac{2}{27}}}{\displaystyle \frac{{\widehat{g}}_3}{8}}-{\displaystyle \frac{1}{4}}\sqrt{a^{-{\displaystyle \frac{1}{27}}}{\displaystyle \frac{{\widehat{g}}_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{a}^{-{\displaystyle \frac{1}{27}}}{\widehat{g}}_2^3}}.$

Now, it is important to note that ${\alpha }_2$ is not another root of the cubic Equation (242)

$0=(\sigma -{\alpha }_1)(\sigma -{\beta }_1)(\sigma -{\gamma }_1)={\left(\sqrt[3]{a}\right)}^2(\sigma -{\alpha }_2)(\sigma -{\beta }_2)(\sigma -{\gamma }_2).$

If one assumes that ${\alpha }_2$ is another root (different from the roots ${\alpha }_1,{\beta }_1,{\gamma }_1$), then this would imply that the polynomial on the left-hand side of (242) should be of the fourth degree, which evidently will be a contradiction. Consequently, the root ${\alpha }_2$ should coincide with one of the roots ${\alpha }_1,{\beta }_1,{\gamma }_1$. Let us first assume that ${\alpha }_2={\alpha }_1$. From the comparison of the corresponding expressions under the cubic and square roots of (235) and (244), it follows that the equality is fulfilled when the following relations between the Weierstrass invariants $g_2,g_3$ and ${\widehat{g}}_2,{\widehat{g}}_3$ are fulfilled:

(245)$g_3=a^{-{\displaystyle \frac{1}{27}}}{\widehat{g}}_3,g_2^3=a^{-{\displaystyle \frac{1}{27}}}{\widehat{g}}_2^3.$

The last relation can be rewritten, expressing the dependance of ${\widehat{g}}_2,{\widehat{g}}_3$ on $g_2,g_3$:

(246)${\widehat{g}}_3=a^{{\displaystyle \frac{1}{27}}}{g}_3,{\widehat{g}}_2=a^{{\displaystyle \frac{1}{81}}}{g}_2.$

It can be noted that these relations are different from the ones obtained in a previous paper [33]:

(247)$g_2^3=a^2{\widehat{g}}_2^3,g_3=-a{\widehat{g}}_3.$

In fact, the correct relations are (245) and (246) (and not those in (247), which are wrong), because they are obtained from the relation for the roots (242). This relation takes into account the multiplication both of the roots and the variables in (242) with the conformal factor $a^{{\displaystyle \frac{2}{9}}}$ and not with $({\sigma }_0,{\sigma }_1,{\sigma }_2)=\sqrt[3]{a}({\sigma }_0,{\sigma }_1,{\sigma }_2)$, as this is applied in the paper [33].

12.2. Finding the Relation Between the Other Two Roots $({\beta }_2,{\gamma }_2)$ and $({\beta }_1,{\gamma }_1)$ for the Case When ${\alpha }_2={\alpha }_1$—The Caseof Positive Discriminant

The other two roots (${\beta }_2,{\gamma }_2)$ of the second Equation (241) $4{\sigma }^3-{\widehat{g}}_2\sigma -{\widehat{g}}_3=0$ can be written similarly to the expressions (239) and (240) for ${\beta }_1$ and ${\gamma }_1$, but with indices “2”:

(248)${\beta }_2=A_2ϵ+B_2{ϵ}^2=-{\displaystyle \frac{(A_2+B_2)}{2}}+i{\displaystyle \frac{\sqrt[.]{3}}{2}}(A_2-B_2),$

(249)${\gamma }_2=A_2{ϵ}^2+B_2{ϵ}^{.}.=-{\displaystyle \frac{(A_2+B_2)}{2}}-i{\displaystyle \frac{\sqrt[.]{3}}{2}}(A_2-B_2).$

Taking into account that ${\alpha }_1={\alpha }_2$, and also the equality (244), it can be obtained that

(250)$A_2+B_2=a^{{\displaystyle \frac{2}{9}}}(A_1+B_1)\Rightarrow A_2-B_2=a^{{\displaystyle \frac{2}{9}}}(A_1-B_1).$

If, substituted in the formulaes (248) and (249) for ${\beta }_2$ and ${\gamma }_2$, the following relations are found between the roots (${\beta }_2,{\gamma }_2)$ and (${\beta }_1,{\gamma }_1)$:

(251)${\beta }_2=a^{{\displaystyle \frac{2}{9}}}{\beta }_1,{\gamma }_2=a^{{\displaystyle \frac{2}{9}}}{\gamma }_1,$

where $({\beta }_1,{\gamma }_1)$ are the complex roots, given by the formulaes (239) and (240).

12.3. Expressing the Ratio of the Weierstrass Invariants $g_2$ and $g_3$ as a Rational Function of Higher-Order Polynomials, Depending on the Modulus Parameter q of the Elliptic Integral in the Legendre Form—The Case of Positive Discriminant

Finding the relation between the Weierstrass invariants $g_2,g_3$ and ${\widehat{g}}_2,{\widehat{g}}_3$ in (245) and (246) is rather insufficient and does not give any valuable information because of the following two reasons:

A. The variable a is not just a parameter, but represents a complicated function, depending on the ratio ${\displaystyle \frac{g_3}{g_2}}$ and on the rational function $K\left(q\right)$, and is determined in (146) as $a={\displaystyle \frac{9g_3}{g_2}}K\left(q\right)$, where (148) $K\left(q\right)\equiv {\displaystyle \frac{q^4-q^2+1}{2q^4-5q^2+2}}$.

B. The invariants ${\widehat{g}}_2,{\widehat{g}}_3$ do not contain any useful information—they appeared after the variable change (233)

${\widehat{g}}_2\equiv {\overline{g}}_2{a}^{{\displaystyle \frac{2}{3}}},{\widehat{g}}_3\equiv {\overline{g}}_{3}a,$

when the equality of the two integrals (234)

$\int {\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-g_2\sigma -g_3}}}={\displaystyle \frac{1}{\sqrt[3]{a}}}\int {\displaystyle \frac{d\sigma }{\sqrt{4{\sigma }^3-{\widehat{g}}_2\sigma -{\widehat{g}}_3}}}.$

was obtained. Although we found the relation between the roots of the polynomials in both sides of (234), it is more important to find the dependance of the invariants $g_2,g_3$ on the invariants ${\overline{g}}_2,{\overline{g}}_3$, expressed through polynomials of the modulus q and calculated in the process of the four-dimensional uniformization, presented in the previous sections.

That is why, now, we shall combine the transformations (245) and (233). We obtain

(252)$g_3=a^{-{\displaystyle \frac{1}{27}}}{\widehat{g}}_3=a^{-{\displaystyle \frac{1}{27}}}{\overline{g}}_{3}a=a^{{\displaystyle \frac{26}{27}}}{\overline{g}}_3,$

(253)$g_2={(a^{-{\displaystyle \frac{1}{27}}})}^{{\displaystyle \frac{1}{3}}}{\widehat{g}}_2=a^{-{\displaystyle \frac{1}{81}}}{\overline{g}}_2{a}^{{\displaystyle \frac{2}{3}}}={\overline{g}}_2{a}^{{\displaystyle \frac{53}{81}}}.$

Taking into account (146) for $a={\displaystyle \frac{9g_3}{g_2}}K\left(q\right)$ with (148) for $K\left(q\right)\equiv {\displaystyle \frac{q^4-q^2+1}{2q^4-5q^2+2}}$ and substituting in the above equalities, it follows that

(254)$g_3={\overline{g}}_3{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{26}{27}}}.{\left({\displaystyle \frac{g_3}{g_2}}\right)}^{{\displaystyle \frac{26}{27}}},$

(255)$g_2={\overline{g}}_2{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{53}{81}}}.{\left({\displaystyle \frac{g_3}{g_2}}\right)}^{{\displaystyle \frac{53}{81}}}.$

Dividing the above two equalities, it can be obtained that

(256)${\displaystyle \frac{g_3}{g_2}}={\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}{\left(9K\left(q\right)\right)}^{\left[{\displaystyle \frac{26.3}{27.3}}-{\displaystyle \frac{53}{81}}\right]}{\left({\displaystyle \frac{g_3}{g_2}}\right)}^{\left[{\displaystyle \frac{26}{27}}-{\displaystyle \frac{53}{81}}\right]}.$

After some transformations, the ratio ${\displaystyle \frac{g_3}{g_2}}$ can be expressed as a function of powers of the ratio ${\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}$ (which is a rational function of polynomials of the modulus parameter q-see (143) and (213)) and powers of $\left(9K\right(q\left)\right)$, which is also a rational function of higher-order polynomials in the nominator and in the denominator according to (148):

(257)${\displaystyle \frac{g_3}{g_2}}={\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}.}.$

This is an important formulae, enabling to express $g_3$ as a linear function of $g_2$ and a complicated function $N\left(q\right)$:

(258)$g_3=g_{2}N\left(q\right),N\left(q\right):={\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}.}.$

In such a way, substituting the above expression in formulae (235) for the (real-valued) first root ${\alpha }_1=A_1+B_1$ of the cubic polynomial $4{\sigma }^3-g_2\sigma -g_3=0$, one can obtain

${\alpha }_1=A_1+B_1=\sqrt[3]{g_2}[\sqrt[3]{F\left(q\right)+{\displaystyle \frac{1}{4}}\sqrt[.]{\tilde{\Delta }}}$

(259)$+\sqrt[3]{F\left(q\right)-{\displaystyle \frac{1}{4}}\sqrt[.]{\tilde{\Delta }}}]=\sqrt[3]{g_2}.\left[{\tilde{A}}_1+{\tilde{B}}_1\right],$

where $F\left(q\right)$ and $\tilde{\Delta }$ (obtained from the discriminant formulae (236) $\Delta ={\displaystyle \frac{g_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{g}_2^3>0$) denote the following functions:

(260)$F\left(q\right):={\displaystyle \frac{1}{8}}{\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}.}={\displaystyle \frac{1}{8}}N\left(q\right),$

(261)$\tilde{\Delta }:={\displaystyle \frac{1}{4}}{\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{28}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{28}}.}-{\displaystyle \frac{1}{4.27}}{g}_2.$

In (259), the following new notations were introduced:

(262)${\tilde{A}}_1:=\sqrt[3]{F\left(q\right)+{\displaystyle \frac{1}{4}}\sqrt[.]{\tilde{\Delta }}},{\tilde{B}}_1:=\sqrt[3]{F\left(q\right)-{\displaystyle \frac{1}{4}}\sqrt[.]{\tilde{\Delta }}}.$

Since the “modified” discriminant $\tilde{\Delta }$ (261) is linearly dependant on the Weierstrass invariant $g_2$, it is evident that $g_2$ can be found if the ratio ${\displaystyle \frac{{\tilde{A}}_1}{{\tilde{B}}_1}}$ is known. This will be used in the next sections.

12.4. A Simple Proof for the Positivity of the Weierstrass Invariants $g_2$ for the One Pair of Weierstrass Invariants ${\overline{g}}_2$ (226) and ${\overline{g}}_3$ (227)

From (261), the condition for a non-zero discriminant can be written as

(263)$g_2\ne 27{\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{28}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{28}}}.$

Now, we shall prove that this expression is different from zero and also is positive. The nominator of the expression for (148) $K\left(q\right)\equiv {\displaystyle \frac{q^4-q^2+1}{2q^4-5q^2+2}}$ and also the expression for ${\overline{g}}_3$ (227),

${\overline{g}}_3=-{\displaystyle \frac{1}{6}}{q}^2(1+q^2)+{\displaystyle \frac{{(1+q^2)}^3}{6^3}},$

are different from zero. For the nominator of $K\left(q\right)$, the positivity is easily established, because it can be represented in the form

(264)$q^4-q^2+1=q^4-2q^2+q^2+1={(q^2-1)}^2+q^2>0,$

which, evidently, is greater than zero.

Consequently, whether $g_2$ will be positive depends on the ratio ${\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}$. The invariant ${\overline{g}}_2$ (226)

${\overline{g}}_2=q^2+{\displaystyle \frac{1}{12}}{(1+q^2)}^2>0$

is evidently positive.

Let us check the sign of expression (227) for ${\overline{g}}_3$, which can be rewritten as

(265)${\overline{g}}_3={\displaystyle \frac{(q^2+1)}{6^3}}(q^4-2{\displaystyle \frac{3.17}{2}}{q}^2+1)$

(266)$={\displaystyle \frac{(q^2+1)}{6^3}}\left[{\left(q^2-{\displaystyle \frac{3.17}{2}}\right)}^2+1-{\left({\displaystyle \frac{3.17}{2}}\right)}^2\right].$

This expression will be positive when the inequalities

(267)$q^2>{\displaystyle \frac{3.17}{2}}+{\displaystyle \frac{\sqrt[.]{6.{17}^2-4}}{2}}>0,$

and also

(268)$q^2>{\displaystyle \frac{3.17}{2}}-{\displaystyle \frac{\sqrt[.]{6.{17}^2-4}}{2}}>0$

are fulfilled. The first case (267) should be disregarded, because it means that $q>6.804$, which is impossible because q should be in the range $0<q<1$.

The second expression (268) is positive when $q>2.168$, which is also impossible. Consequently, the expression in the square brackets in (266) cannot be positive when $q>6.804$ and $q>2.168$.

However, when $q<6.804$ and $q<2.168$ (always fulfilled since $0<q<1$), the expression

(269)$\left(q^2-{\displaystyle \frac{3.17}{2}}-{\displaystyle \frac{\sqrt[.]{6.{17}^2-4}}{2}}\right)\left(q^2-{\displaystyle \frac{3.17}{2}}+{\displaystyle \frac{\sqrt[.]{6.{17}^2-4}}{2}}\right)>0$

is positive, because both expressions in the round brackets are negative and (269) is with a plus sign.

Therefore, we proved that for all values of q, the invariant $g_2$ (263) is positive and from the assumption for the positivity of the discriminant (236)

$\Delta ={\displaystyle \frac{g_3^2}{4}}-{\displaystyle \frac{1}{108}}{g}_2^3>0\Rightarrow g_3^2>{\displaystyle \frac{1}{27}}{g}_2^3,$

it follows also that $\mid g_3\mid >{\displaystyle \frac{1}{3\sqrt{3}}}\mid g_2^{{\displaystyle \frac{3}{2}}}\mid$.

12.5. The Possibility for a Positive and Real Invariant $g_2$,Imaginary Invariant $g_3$ and a Negative Discriminant—A Theorem in[60]

Formally, it may seem that if $g_3$ is negative but the condition for the moduli is fulfilled, then the discriminant $\Delta$ will be positive. In fact, here, we may have the case when $g_2$ is positive but $g_3$ can be imaginary. This problem is not investigated in the mathematical literature, but it is important that the possibility for $g_3$ to be imaginary does not contradict a theorem in problem 12, paragraph 6 of chapter 1 in the monograph by Koblitz [60]:

Theorem 12. 

The roots $e_1,e_2,e_3$ of the polynomial $4x^3-g_{2}x-g_3=0$ are all real ones if and only if the invariants $g_2$ and $g_3$ are all real ones and the discriminant $\Delta =g_2^3-27g_3^2$$>0$ is positive.

In this paper, we are not investigating the case when all the roots of the polynomial $4x^3-g_{2}x-g_3=0$ are real; we deal with the case when only the first root ${\alpha }_1$ is a real one and the other two roots ${\beta }_1$ and ${\gamma }_1$ are imaginary. Consequently, it is not obligatory that the invariants $g_2$ and $g_3$ should be real ones. We obtained that the invariant $g_2$ is positive, but let us remember that for the invariant $g_3$, we found the formulae (258)

$g_3=g_{2}N\left(q\right),N\left(q\right):={\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}}.$

We already proved that the ratio ${\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}$ is positive and real; also, we proved that the nominator of the function (148) $K\left(q\right)\equiv {\displaystyle \frac{q^4-q^2+1}{2q^4-5q^2+2}}$ is positive. The denominator can be decomposed as

(270)$2q^4-5q^2+2=2\left[\left(q^2-2\right)\left(q^2-{\displaystyle \frac{1}{2}}\right)\right].$

Now, if one assumes that this expression is positive, then this is possible only when the following inequalities are simultaneously fulfilled:

(271)$q<\sqrt{2}\approx 1.41,q<{\displaystyle \frac{1}{\sqrt{2}}}\approx 0.707106.$

The other two inequalities with the positive sign should be discarded because $0<q<1$. But if the inequalities (271) are fulfilled and $K\left(q\right)$, $N\left(q\right)$ and $g_3$ are positive, then from the real-valuedness of $g_2$ and $g_2$ and also the positivity and real-valuedness of the discriminant $\Delta =g_2^3-27g_3^2$$>0$ from the theorem in [60], it should follow that the roots $e_1,e_2,e_3$ of the polynomial $4x^3-g_{2}x-g_3=0$ should all be real. But this will be in contradiction to the case of one real root and two imaginary roots, investigated in this paper.

Consequently, an important conclusion can be made that there can be two cases:

1st case: positive invariant $g_2$, imaginary $g_3$ (or both imaginary invariants $g_2$ and $g_3$) and positive discriminant $\Delta =g_2^3-27g_3^2$$>0$. Then, not all the roots should be real ones (the presently investigated case).

2nd case: real roots ${\alpha }_1,{\beta }_1,{\gamma }_1$ (the case not investigated in this paper). Then, the invariants $g_2$ and $g_3$ should be real and the discriminant $\Delta$ should be negative. From the point of view of practical applications for the inter-satellite communications, it is evident that the case of a negative discriminant and modulus parameter $q<{\displaystyle \frac{1}{\sqrt{2}}}\approx 0.707106$ (remember that for a GPS orbit $e\approx 0.01$ and $q={\displaystyle \frac{1-e}{1+e}}$ ), such a case seems to be realistic and should be investigated.

Note that in the strict mathematical sense, the theorem in the monograph by Koblitz [60] should be reformulated and proved.

The case of a negative discriminant for a general cubic equation $x^3+px+\tilde{q}=0$ is investigated in details in the monograph by Obreshkoff [61] and the roots have been expressed by the concise formulae

(272)$x_{0,1,2}=2\sqrt{-{\displaystyle \frac{p}{3}}}cos{\displaystyle \frac{\phi +2k\pi }{3}},k=0,1,2$

and the angular variable $\phi$ is given by the formulae

(273)$tg\phi =-{\displaystyle \frac{2\sqrt{-\Delta }}{\tilde{q}}}.$

We have not investigated the case of the other pair of invariants (212) ${\overline{g}}_2^{(1,2)}$ (identical with the invariant (226)) and (214) ${\overline{g}}_3^{(1,2)}$; when the second invariant is given by a polynomial of the sixth degree and depending on its roots, it can be either positive or negative.

12.6. Finding the Weierstrass Function $\rho \left(z\right)$ and theWeierstrass Invariant $g_2$ From the Known Roots of the Cubic Equation andthe Group-Theoretical Law for Summing up of Points on the Cubic Curve

The main goal of every investigation, related to the cubic algebraic equation in its parametrizable form, is to find the expression for the Weierstrass function $\rho \left(z\right)$ by known Weierstrass invariants $g_2$ and $g_3$. Since we started from a Weierstrass integral, derived from an elliptic integral in the Legendre form, previously, we found how the invariant $g_3$ can be expressed through the other invariant $g_2$ (258)

$g_3=g_{2}N\left(q\right),N\left(q\right):={\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}.}.$

Now, it remains to find the other Weierstrass invariant $g_2$ and the Weierstrass function itself. This shall be performed by using the found roots ${\alpha }_1,{\beta }_1,{\gamma }_1$ for the case of a positive discriminant (236) $\Delta ={\displaystyle \frac{g_3^2}{4}}-{\displaystyle \frac{1}{27.4}}{g}_2^3>0$, the relation (274) between these roots with the values of the Weierstrass function at the points on the complex plane ${\omega }_1,{\omega }_2,{\omega }_3={\omega }_1+{\omega }_2$ (for convenience we omit the factor ${\displaystyle \frac{1}{2}}$ in the argument of the Weierstrass function $\rho \left(z\right)$),

(274)$e_1=\rho \left({\omega }_1\right)={\alpha }_1,e_2=\rho \left({\omega }_2\right)={\beta }_1,e_3=\rho ({\omega }_1+{\omega }_2)={\gamma }_1,$

and also the group-theoretical law for summation of the points ${\omega }_1,{\omega }_2,{\omega }_3={\omega }_1+{\omega }_2$ on the elliptic curve, expressed analytically by the formulae (199)

$\rho (z_1+z_2)=\rho \left(z_1\right)+\rho \left(z_2\right)+{\displaystyle \frac{1}{4}}{\left({\displaystyle \frac{{\rho }^{\prime }\left(z_1\right)-{\rho }^{\prime }\left(z_2\right)}{\rho \left(z_1\right)-\rho \left(z_2\right)}}\right)}^2.$

Taking into account (274), let us take the ratio ${\displaystyle \frac{\rho \left(z_1\right)}{\rho \left(z_2\right)}}$ of the Weierstrass function at the two different values $z_1={\omega }_1$ and $z_2={\omega }_2$ and also the ratio ${\displaystyle \frac{\rho \left(z_1\right)}{\rho \left(z_3\right)}}$ at the two values $z_1={\omega }_1$ and $z_3={\omega }_3={\omega }_1+{\omega }_2$:

(275)${\displaystyle \frac{\rho \left(z_1\right)}{\rho \left(z_2\right)}}={\displaystyle \frac{{\alpha }_1}{{\beta }_1}},{\displaystyle \frac{\rho \left(z_1\right)}{\rho \left(z_3\right)}}={\displaystyle \frac{{\alpha }_1}{{\gamma }_1}},$

where the root ${\alpha }_1=A_1+B_1$ is given by formulae (235) and the complex roots ${\beta }_1$ and ${\gamma }_1$ are given by formulaes (239) and (240), respectively. Since ${\beta }_1$ is a complex root, the Weierstrass function $\rho \left(z_2\right)={\beta }_1$ is also a complex function and can be represented as $\rho \left(z_2\right)={\rho }_1\left(z_2\right)+i_2\rho \left(z_2\right)$, where ${\rho }_1\left(z_2\right)$ is the real part of $\rho \left(z_2\right)$ and ${\rho }_2\left(z_2\right)$ is the imaginary part.

12.6.1. The First Equation

The first equation in (275) can be written as

(276)${\displaystyle \frac{\rho \left(z_1\right)}{{\rho }_1\left(z_2\right)+i{\rho }_2\left(z_2\right)}}={\displaystyle \frac{A_1+B_1}{\frac{-(A_1+B_1)}{2}+i\frac{\sqrt[.]{3}}{2}(A_1-B_1)}}.$

Denoting

(277)$m:={\displaystyle \frac{{\tilde{A}}_1-{\tilde{B}}_1}{{\tilde{A}}_1+{\tilde{B}}_1}},$

the above equation can be represented as

(278)$-{\displaystyle \frac{{\rho }_{.}\left(z_1\right)}{2}}-{\rho }_1\left(z_2\right)+i{\displaystyle \frac{\sqrt[.]{3}}{2}}{\rho }_{.}\left(z_1\right)m-i{\rho }_2\left(z_2\right)=0.$

Setting up equal to zero the real and imaginary part, one can obtain

(279)${\rho }_{1.}\left(z_2\right)=-{\displaystyle \frac{{\rho }_{.}\left(z_1\right)}{2}},{\rho }_{2.}\left(z_2\right)={\displaystyle \frac{\sqrt[.]{3}}{2}}{\rho }_{.}\left(z_1\right)m.$

12.6.2. The Second Equation

Let us now make use of the second equation in (275), which can be written as

(280)${\displaystyle \frac{\rho \left(z_1\right)}{{\rho }_{.}\left(z_3\right)}}={\displaystyle \frac{A_1+B_1}{-\frac{(A_1+B_1)}{2}-i\frac{\sqrt[.]{3}}{2}(A_1-B_1)}}.$

Using the group-theoretical law for summing up of points, represented by formulae (199) for expressing ${\rho }_{.}\left(z_3\right)={\rho }_{.}(z_1+z_2)$, the above equation can be written as

${\displaystyle \frac{\rho \left(z_1\right)}{\rho \left(z_1\right)+\rho \left(z_2\right)+\frac{1}{4}{\left(\frac{{\rho }^{\prime }\left(z_1\right)-{\rho }^{\prime }\left(z_2\right)}{\rho \left(z_1\right)-\rho \left(z_2\right)}\right)}^2}}$

(281)$={\displaystyle \frac{1}{-\frac{1}{2}-i\frac{\sqrt[.]{3}}{2}m}}.$

Using the two equalities in (279), Equation (281) can be represented as $Q_1+i{Q}_2=0$, where the real part $Q_1$ and the complex part $Q_2$ are

$Q_1:=-3{\rho }^3\left(z_1\right)(5m^2+3)-{\displaystyle \frac{9}{4}}{\left({\rho }^{\prime }\left(z_1\right)\right)}^2(1-{\displaystyle \frac{1}{3}}{m}^2)$

(282)$+{\displaystyle \frac{3}{4}}\rho \left(z_1\right){\displaystyle \frac{\partial m}{\partial q}}.{\displaystyle \frac{\partial q}{\partial z}}+{\displaystyle \frac{3}{2}}\rho \left(z_1\right){\rho }^{\prime }\left(z_1\right)m{\displaystyle \frac{\partial m}{\partial q}}.{\displaystyle \frac{\partial q}{\partial z}}=0,$

(283)$Q_2:={\rho }^3\left(z_1\right)m(m^2-1)+{\displaystyle \frac{1}{2}}{\rho }^{\prime }\left(z_1\right)\left[{\rho }^{\prime }\left(z_1\right)m+\rho \left(z_1\right){\displaystyle \frac{\partial m}{\partial q}}.{\displaystyle \frac{\partial q}{\partial z}}\right]=0.$

Since the equation $Q_1+i{Q}_2=0$ is fulfilled when the real and imaginary part are equal to zero, we have to set up $Q_1=0$ and $Q_2=0$. Also, we have assumed that the modulus parameter q of the elliptic integral in the Legendre form might depend on the complex coordinate z, on which, by definition, depends the Weierstrass function $\rho \left(z\right)$. However, in view of the complexity of the equations, for now, we shall consider the simplified case ${\displaystyle \frac{\partial q}{\partial z}}=0$.

12.6.3. First Way for Representation of $m^2$ and of ${\left({\rho }^{\prime }\left(z_1\right)\right)}^2$

Further, from the second Equation (283) for $Q_2$, one may express the ratio (277) $m^2:={\left({\displaystyle \frac{{\tilde{A}}_1-{\tilde{B}}_1}{{\tilde{A}}_1+{\tilde{B}}_1}}\right)}^2$ as

(284)$m^2={\displaystyle \frac{{\rho }^3\left(z_1\right)-\frac{1}{2}{\left({\rho }^{\prime }\left(z_1\right)\right)}^2}{{\rho }^3\left(z_1\right)}}.$

Substituting $m^2$ from (284) in the first Equation (282) for $Q_1$ and denoting ${\left({\rho }^{\prime }\left(z_1\right)\right)}^2:=X$, one obtains the following quadratic equation with respect to the variable X:

(285)$9X^2-18.8{\rho }^3\left(z_1\right)X+72{\rho }^6\left(z_1\right)=0.$

It is known that the roots of the simple quadratic equation $x^2+px+q=0$ are found from the formulae [61]

(286)$x=-{\displaystyle \frac{p}{2}}\pm \sqrt[.]{{\displaystyle \frac{p^2}{4}}-q}.$

Consequently, the roots of the quadratic Equation (285) are

(287)${\left({\rho }^{\prime }\left(z_1\right)\right)}^2:=X=18.4{\rho }^3\left(z_1\right)\pm \sqrt[.]{{\displaystyle \frac{{\left(18.8\right)}^2{\rho }^6\left(z_1\right)}{4}}-72}.$

This relation can be substituted in formulae (284) for $m^2$ and for the moment investigating only the case of a positive sign in front of the square root, $m^2$ can be written as

(288)$m^2=-35-3\sqrt[.]{{12}^2-{\displaystyle \frac{2}{{\rho }^6\left(z_1\right)}}}.$

Respectively, the solution (287) for ${\left({\rho }^{\prime }\left(z_1\right)\right)}^2$ can be represented as

(289)${\left({\rho }^{\prime }\left(z_1\right)\right)}^2=72{\rho }^3\left(z_1\right)+6{\rho }^3\left(z_1\right)\sqrt[.]{{12}^2-{\displaystyle \frac{2}{{\rho }^6\left(z_1\right)}}}.$

The last two relations are the first way for representing $m^2$ and ${\left({\rho }^{\prime }\left(z_1\right)\right)}^2$. Now, it shall be proved that there is a second way for representing these expressions.

12.6.4. Second Way for Representation of $m^2$ and of ${\left({\rho }^{\prime }\left(z_1\right)\right)}^2$

The second way is based on the standard formulae for the parametrizable form of a cubic curve:

(290)${\left({\rho }^{\prime }\left(z_1\right)\right)}^2=4{\rho }^3\left(z_1\right)-g_2\rho \left(z_1\right)-g_3.$

Substituting into (284) for $m^2$, we obtain the second representation:

(291)$m^2=-2+{\displaystyle \frac{g_2}{2}}.{\displaystyle \frac{1}{{\rho }^2\left(z_1\right)}}+{\displaystyle \frac{g_3}{2}}{\displaystyle \frac{1}{{\rho }^3\left(z_1\right)}}.$

Next we shall use also expression (258) for $g_3$ through $g_2$:

$g_3=g_{2}N\left(q\right),N\left(q\right):={\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}}.$

Taking this equality into account and also that according to (260), $N\left(q\right)=8F\left(q\right)$, the second representation (291) for $m^2$ can be written also in the form

(292)$m^2={\displaystyle \frac{-{\rho }^3\left(z_1\right)+\left[\frac{\rho \left(z_1\right)}{2}+4F\left(q\right)\right]g_2}{{\rho }^3\left(z_1\right)}}.$

Setting up equal the two representations (288) and (291) for $m^2$, one can obtain

$-34{\rho }^3\left(z_1\right)-3{\rho }^3\left(z_1\right)\sqrt[.]{{12}^2-{\displaystyle \frac{2}{{\rho }^6\left(z_1\right)}}}$

(293)$={\displaystyle \frac{g_2}{2}}(\rho \left(z_1\right)+8F\left(q\right)).$

12.6.5. Another Interpretation of the Equation (293)for the Equality of the Two Representations of $m^2$—A Theorem fromHigher Algebra for Investigating Whether There are Roots of an AlgebraicPolynomial

After some transformations, the above equation can be represented in the form of a sixth-order algebraic equation with respect to $\rho \left(z_1\right)$:

$\left[9.{12}^2\right]{\rho }^6\left(z_1\right)-34g_2{\rho }^4\left(z_1\right)-68.4F\left(q\right){\rho }^3\left(z_1\right)$

(294)$-g_2^2\left[{\rho }^2\left(z_1\right)-4F\left(q\right)\rho \left(z_1\right)+16F^2\left(q\right)\right]=0.$

This equation of the 6th degree is interesting, because it does not always have roots. It will, however, have roots for certain values of the invariant $g_2$. Thus, the values of $\rho \left(z_1\right)$ are interrelated to the values of $g_2$. It is interesting to investigate the following problem: According to (293), for a given function $\rho \left(z_1\right)$, there should be only one value of $g_2$. However, the algebraic Equation (294), if solved as a quadratic equation with respect to $g_2$ and for a given function $\rho \left(z_1\right)$, should have two roots with respect to $g_2$, which will contradict the unique expression of $g_2$ from (293).

Now, let us give the formulation of the theorem of Schur from higher algebra, which gives the possibility to establish whether a polynomial has roots. The formulation of the theorem, originally published by Schur [62] in 1918, is taken from the monograph [61], where many other theorems concerning roots of polynomials can be found. The theorem of Schur has been applied in the paper [20] for establishing whether two different algebraic equations of the fourth degree have roots.

The first of these equations was related to the introduced in [20] physical notion of a “space-time interval”, which, for a particular case, was proved to be positive, negative or even equal to zero. By means of the Schur theorem, it was proved that the corresponding algebraic equation has roots.

The second equation was related to the notion of “geodesic distance”, which is related to the propagation of a light or radio signal in a gravitational field and, therefore, should be positive. Respectively, again, by applying the Schur theorem, it was proved that the corresponding algebraic equation does not have roots. The confirmation of the physical definitions has been commented on in the paper [40] and also the review paper [22]. This consistency between the mathematical results and the corresponding physical definitions give some grounds to believe that the Schur theorem is a reliable mathematical tool for the investigation, also, of Equation (294).

Theorem 13 

([62]). The necessary and sufficient conditions for the polynomial of n-th degree

(295)$f\left(y\right)=a_0{y}^n+a_1{y}^{n-1}+....+a_{n-2}{y}^2+a_{n-1}y+a_n$

to have roots only in the circle $\mid y\mid <1$ are the following ones: 1. 

The fulfillment of the inequality

(296)$\mid a_0\mid >\mid a_n\mid .$

2. 

The roots of the polynomial of $(n-1)$ degree

(297)$f_1\left(y\right)={\displaystyle \frac{1}{y}}\left[a_{0}f\left(y\right)-a_n{f}^{*}\left(y\right)\right]$

should be contained in the circle $\mid y\mid <1$, where $f^{*}\left(y\right)$ is the s.c. “inverse polynomial”, defined as

(298)$f^{*}\left(y\right)=y^{n}f\left({\displaystyle \frac{1}{y}}\right)=a_n{y}^n+a_{n-1}{y}^{n-1}+....+a_2{y}^2+a_{1}y+a_0.$

In the case of fulfillment of the inverse inequality

(299)$\mid a_0\mid <\mid a_n\mid ,$

the $(n-1)$ degree polynomial $f_1\left(y\right)$ (again with the requirement the roots to remain within the circle $\mid y\mid <1$) is given by the expression

(300)$f_1\left(y\right)=a_{n}f\left(y\right)-a_0{f}^{*}\left(y\right).$

In the paper [20], the Schur theorem has been generalized for a “chain” of algebraic equations, when it is applied to a chain of algebraic equations with diminishing degrees. Many other theorems for investigating when a polynomial has roots or does not have can be found also in the contemporary monograph [69].

12.7. Finding the Weierstrass Function $\rho \left(z_1\right)$ and theInvariant $g_2$

Equation (293), further, will play an important role, since it allows expressing the invariant $g_2$ as a function of the Weierstrass function $\rho \left(z_1\right)$. More importantly, it shows that for a given function $\rho \left(z_1\right)$, there is only one value of $g_2$.

Yet, Equation (293) is not sufficient for finding both the Weierstrass invariant $g_2$ and the Weierstrass function $\rho \left(z_1\right)$. We need one more equation, and for the purpose, we will use the defining equality for m (277) $m:={\displaystyle \frac{{\tilde{A}}_1-{\tilde{B}}_1}{{\tilde{A}}_1+{\tilde{B}}_1}}$, where ${\tilde{A}}_1$ and ${\tilde{B}}_1$ were previously defined by the formulaes (262)

${\tilde{A}}_1:=\sqrt[3]{F\left(q\right)+{\displaystyle \frac{1}{4}}\sqrt[.]{\tilde{\Delta }}},{\tilde{B}}_1:=\sqrt[3]{F\left(q\right)-{\displaystyle \frac{1}{4}}\sqrt[.]{\tilde{\Delta }}}.$

Thus, it can easily be obtained that ${\displaystyle \frac{{\tilde{A}}_1}{{\tilde{B}}_1}}={\displaystyle \frac{1+m}{1-m}}$. Taking into account (262), it can be derived that

(301)${\displaystyle \frac{F\left(q\right)+\frac{1}{4}\sqrt[.]{\tilde{\Delta }}}{F\left(q\right)-\frac{1}{4}\sqrt[.]{\tilde{\Delta }}}}={\displaystyle \frac{{(1+m)}^3}{{(1-m)}^3}},$

where $\tilde{\Delta }$ has been defined by the equality (261)

$\tilde{\Delta }:={\displaystyle \frac{1}{4}}{\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{28}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{28}}}-{\displaystyle \frac{1}{4.27}}{g}_2.$

Denoting the first term in the above equality as ${\tilde{\Delta }}_1:=$ ${\displaystyle \frac{1}{4}}{\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{28}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{28}}}$ and remembering the defining equality (260)

$F\left(q\right):={\displaystyle \frac{1}{8}}{\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}}={\displaystyle \frac{1}{8}}N\left(q\right),$

one can obtain the simple expression

(302)${\displaystyle \frac{{\tilde{\Delta }}_1}{F\left(q\right)}}=2{\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}.}=16F\left(q\right).$

Consequently, the discriminant $\tilde{\Delta }$ (261) can be represented as

(303)$\tilde{\Delta }:={\tilde{\Delta }}_1-{\displaystyle \frac{1}{4.27}}{g}_2=16F^2\left(q\right)-{\displaystyle \frac{1}{4.27}}{g}_2.$

Making use of (301), one can find an expression for $\sqrt[.]{\tilde{\Delta }}$:

(304)$\sqrt[.]{\tilde{\Delta }}=4F\left(q\right){\displaystyle \frac{\left[{(1+m)}^3-{(1-m)}^3\right]}{\left[{(1+m)}^3+{(1-m)}^3\right]}}.$

Transforming this equality and using (303), one obtains the final formulae for $\tilde{\Delta }$, from where $g_2$ can be expressed and $m^2$ can be expressed in two different ways:

(305)$\tilde{\Delta }=16F^2\left(q\right)-{\displaystyle \frac{1}{4.27}}{g}_2={\displaystyle \frac{16F^2\left(q\right)m^2{(3+m^2)}^2}{{(1+3m^2)}^2}}.$

In this formulae, we substitute $m^2$, expressed in (291) through the Weierstrass invariant $g_2$ and higher powers of the Weierstrass function $\rho \left(z\right)$. After some lengthy calculations, the following complicated expression is obtained:

(306)${\displaystyle \frac{1}{4.27}}{g}_2=16F\left(q\right)+{\displaystyle \frac{16F^2\left(q\right){\rho }^6\left(z_1\right)}{{\rho }^9\left(z_1\right)}}{\displaystyle \frac{P_1(\rho \left(z_1\right),q)}{P_2(\rho \left(z_1\right),q)}},$

where $P_1(\rho \left(z_1\right),q)$ and $P_2(\rho \left(z_1\right),q)$ are the expressions (note that, everywhere, the dot “.” means multiplication, for example, $35.{32}^2=35\times {32}^2$)

$P_1(\rho \left(z_1\right),q):=-35.{32}^2{\rho }^9\left(z_1\right)+32.306{\rho }^6\left(z_1\right)\sqrt{{12}^2{\rho }^6\left(z_1\right)-2}$

(307)$+\left(99.9\right){\rho }^3\left(z_1\right)\left({12}^2{\rho }^6\left(z_1\right)-2\right)+27{({12}^2{\rho }^6\left(z_1\right)-2)}^{{\displaystyle \frac{3}{2}}},$

$P_2(\rho \left(z_1\right),q):={104}^2{\rho }^6\left(z_1\right)+\left(104.18\right){\rho }^3\left(z_1\right)\sqrt{12{\rho }^6\left(z_1\right)-2}$

(308)$+81(12{\rho }^6\left(z_1\right)-2).$

This equation and also Equation (293)

$-34{\rho }^3\left(z_1\right)-3{\rho }^3\left(z_1\right)\sqrt[.]{{12}^2-{\displaystyle \frac{2}{{\rho }^6\left(z_1\right)}}}$

$={\displaystyle \frac{g_2}{2}}(\rho \left(z_1\right)+8F\left(q\right))$

can be considered as a system of two equations with respect to $g_2$ and $\rho \left(z_1\right)$. If, from (293), the Weierstrass invariant $g_2$ is expressed and then substituted into (306), then a complicated equation is obtained only with respect to $\rho \left(z_1\right)$. It is not clear how to find its roots (and also how many there are), but if $\rho \left(z_1\right)$ is found, then $g_2$ is uniquely expressed from (293).

There is one more possibility, namely, to use the second representation for $m^2$ by formulae (292), which is to be substituted in Equation (305) for $\tilde{\Delta }$. We shall not present here this method, which is more complicated and will lead to the derivation of a cubic equation with respect to $g_2$, which, after solving, has to be substituted in the earlier derived algebraic Equation (294) of the sixth order with respect to $\rho \left(z_1\right)$. However, again, a higher-order algebraic equation has to be solved, so this method is not effective.

13. Conclusions

In this paper, the propagation time of a signal was analytically calculated by making use of the null cone equation $d{s}^2=0$ from the metric element, defined for the near-Earth space. Two different cases were considered, corresponding to two different parametrizations of the space coordinates $(x,y,z)$ of the satellite orbit: The first parametrization is in terms of the eccentric anomaly angle E, which is the dynamical parameter for the case of a plane elliptic orbit, and the second parametrization is in terms of all the six Kepler parameters $(f,a,e,\Omega ,I,\omega )$ for the case of space-oriented orbit with only the true anomaly angle f considered to be the dynamical parameter for the motion of the satellite.

The propagation time for the first case is proved to contain the three known elliptic integrals of the first, second and the third kind; all of them are zero-order elliptic integrals. The definitions of the different kinds of elliptic integrals are given in Section 4.3.

The propagation time for the second case is investigated in detail in Section 6. It is expressed in more complicated elliptic integrals of the second and of the fourth order. These are higher-order elliptic integrals, the mathematical definitions of which are given in Section 4.1. It should be mentioned these integrals have been known in the mathematical literature for a long time, but the known physical applications are not so many compared to the elliptic integrals of the first, second and third kind.

For both cases, the sum of the $O\left(c^{-1}\right)$ and $O\left(c^{-3}\right)$ propagation time corrections are real-valued expressions. For the case of plane elliptical orbits, this was proved in Section 5.1, and for the second case of space-oriented orbits, this was proved in Section 7. This proves the correctness of the theoretical approach of parametrizing the space coordinates $x,y$ or $x,y,z$ in the null cone equation with the orbital coordinates, known from celestial mechanics.

Numerical calculations for the propagation time for the first case of plane elliptical orbit are performed in Section 8 and they confirm that during the propagation time, the signal travels a much greater distance than the distance travelled by the satellite. It turned out also that theorems from general relativity theory [19] and from differential geometry [24,26] confirm the conclusion that various parametrizations of the space coordinates, related to the orbital motion of the satellite, can be used. This proof comes in three stages: The real-valuedness of all the expressions for the propagation time for all the two cases, the theorems from general relativity theory and differential geometry and, lastly, the numerical calculations confirm the substantial fact that no matter whether a parametrization of the space coordinates (related to the satellite orbits) is used, the propagation time for a signal, expressed from the null cone equation, gives the typical values for signal propagation and not for the moving satellite.

It is important to clarify that the calculated propagation time gives the correspondence between the initial and final values of the dynamical parameter (the eccentric anomaly angle E or the true anomaly angle f) and the initial $T_{init}$ and final $T_{fin}$ moments of propagation of the signal. However, in such a theoretical setup, $T_{init}$ and $T_{fin}$ are not related to any real processes of transmission of signals between satellites on one and the same orbit or on different space-oriented orbits. In such a modified setting, the initial position of the signal-emitting satellite and the final position of the signal-receiving satellite have to be given and they will correspond to the initial $T_{init}$ and final $T_{fin}$ moments of propagation of the signal. They have to be calculated in view of the fact that the trajectory of the signal is curved due to the action of the gravitational field. Due to this, the propagation time will not be the (presumed to be known at each moment of time) Euclidean distance between the satellites, divided by the velocity of light. The nontrivial moment in such a problem is that both satellites are moving (at the moments of emission and receiving of the signal and also during the time of propagation of the signal) so that the moments of emission and reception of the signal (they are denoted by $T_1$ and $T_2$) should incorporate these important peculiarities. In other words, the curvature of the signal trajectory will have to be such that at a certain space point of the second satellite, the signal will be intercepted, provided also that the signal is emitted at the proper position by the first (signal emitting) satellite.

The last three sections of this paper (Section 10.1, Section 11.1 and Section 12) are dedicated to the problem of finding new solutions of elliptic integrals in the Legendre form by means of transforming them to elliptic integrals in the Weierstrass form, which, after inversion, have as a solution the Weierstrass elliptic function. The necessity of studying such integrals (at first, in the pure mathematical aspect, but subsequently, applications in problems related to physics might be searched out) is dictated not only by their important application in GPS inter-satellite communications with account also of the general relativity effects, but also by the wide application of elliptic integrals in various problems of mechanics, cosmology, black hole physics, integrable systems, etc. A short review of these different applications has been given in the paper [33].

In this paper (as well in [33]), the representation of the integral is exposed in two different ways; one of them is based on the transformation $x={\displaystyle \frac{a}{y^2}}+b$ ($a,b$ are proved to be functions of the modulus parameter q of the integral in the Legendre form), and the second method is related to the s.c. “four-dimensional uniformization”, exposed in the monograph by Whittaker and Watson [35]. After comparison of the Weierstrass invariants in the two representations, the main result is the representation (257) for the ratio of the two invariants

${\displaystyle \frac{g_3}{g_2}}={\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right)}^{{\displaystyle \frac{81}{56}}}{\left(9K\left(q\right)\right)}^{{\displaystyle \frac{25}{56}}},$

where $K\left(q\right)$ is the function (148), $K\left(q\right)\equiv {\displaystyle \frac{q^4-q^2+1}{2q^4-5q^2+2}}$, ${\overline{g}}_2$ and ${\overline{g}}_3$ are the Weierstrass invariants, obtained after the four-dimensional uniformization [35]. This result is not trivial and previously not investigated in the monographs on elliptic functions. Further, it was proved also that the Weierstrass invariant $g_2$ and the Weierstrass function $\rho \left(z\right)$ can be found from a complicated system of two nonlinear equations.

Another interesting result was obtained in this paper, suggested by a theorem in the monograph by Koblitz [60]. If $g_2$ is a real number (or a real and positive function) and there is one real root and two imaginary ones for the cubic polynomial $4x^3-g_{2}x-g_3$ (the case investigated in this paper), then the other invariant $g_3$ should be imaginary and the discriminant $\Delta =g_2^3-27g_3^2$ should be negative (a case not investigated in this paper). It should be kept in mind also that the case “negative discriminant and imaginary $g_3$” does not follow directly from the formulation of the theorem, but rather than from the “non-fulfillment” of the theorem, which is an unusual moment. No doubt, more strict mathematical treatment is necessary.

The case of a negative discriminant would correspond to three real roots of the cubic polynomial (see [61,68]) and this case needs a separate investigation. However, the theorem in [60] has a very specific and intricate formulation: From it, for example, follows that for a positive discriminant, one may have also a case with three real roots of the polynomial (contrary to the standard case in the algebraic books [61,68]), but then, the invariants $g_2$ and $g_3$ should also be real! So, the interesting fact is that the Weierstrass invariants $g_2$ and $g_3$ have a purely “algebraic” meaning whether the roots are only real or one real and two complex. Consequently, the theorem in [60] should be reformulated for all the possible cases and investigated again.

Transformations of the elliptic integral with a general fourth-order polynomial under the root in the denominator into an integral in the Legendre form have been studied also in the monograph [70] (also in [42]), but in [70], transformation to an integral in the Weierstrass form has not been performed. However, in [42], the inverse transformation from an integral in the Weierstrass form with the polynomial

(309)$Y^2=(x-\alpha )(x-\beta )(x-\gamma )$

to the integral in the Legendre form [42]

(310)$\sqrt{\alpha -\gamma }\int {\displaystyle \frac{dx}{Y}}=\int {\displaystyle \frac{2dy}{\sqrt{(1-y^2)(1-{\tilde{q}}^2{y}^2)}}},$

was made, where the transformation

(311)$x-\alpha ={\displaystyle \frac{(\alpha -\beta )y^2}{(1-y^2)}}$

has been applied and the modulus of the elliptic integral $\tilde{q}$ in the Legendre form was determined by the formulae

(312)${\tilde{q}}^2={\displaystyle \frac{\beta -\gamma }{\alpha -\gamma }}.$

An interesting problem for future research is the following: The ratio (257) of the invariants $g_2$ and $g_3$ can be used for calculating the roots $\alpha ,\beta ,\gamma$. The resulting expressions will depend on the invariant $g_2$ and can be substituted in the formulae (312) for ${\tilde{q}}^2$. Then, it is interesting to see whether the invariant $g_2$ can be determined in such a way that the calculated ${\tilde{q}}^2$ modulus according to (312) will coincide with the initially given q in the integral $\int {\displaystyle \frac{dy}{\sqrt{(1-y^2)(1-q^2{y}^2)}}}$.

There are also another problems, related to the implementation of the Schur theorem from higher algebra. Besides the ratio (257), another ratio

(313)${\displaystyle \frac{g_3}{g_2}}=\left({\displaystyle \frac{{\overline{g}}_3}{{\overline{g}}_2}}\right){\left(9K\left(q\right)\right)}^{-{\displaystyle \frac{1}{3}}.}$

can be obtained, if it is assumed that the roots ${\alpha }_1$ and ${\alpha }_2$ do not coincide ( ${\alpha }_1={\alpha }_2$), but the roots ${\beta }_1$ and ${\beta }_2$ coincide, i.e., ${\beta }_1={\beta }_2$ This case also is not considered in this paper, but the investigation should be performed in the same manner.

Also, an important problem for further investigation is the algebraic treatment of the higher-order algebraic equation for the discriminant (236) $\Delta ={\displaystyle \frac{g_3^2}{4}}-{\displaystyle \frac{1}{108}}{g}_2^3$ of the cubic polynomial in the Weierstrass elliptic integral, which, in this investigation, was assumed to be positive.

Data Availability Statement

No new data were created or analyzed in this study. Data sharing is not applicable to this article.

Acknowledgments

The author is grateful also to all the three referees of this paper for their constructive remarks. The author is also grateful to the organizers of the Sixteenth Conference of the Euro-American Consortium for Promoting the Application of Mathematics in Technical and Natural Sciences (AMITaNS24) 21–26 June 2024, Albena resourt, Bulgaria (also the previous AMITaNS-conferences in 2022 and 2023), and especially to Mihail Todorov (Technical University, Sofia, Bulgaria; Institute of Mathematics and Informatics, Bulgarian Academy of Sciences, Sofia, Bulgaria) for the opportunity to participate in these conferences and to present some of the material in this paper. The author is grateful also to Sergei Kopeikin (University of Missouri-Columbia, USA) for sparking my interest in Applied General Relativity and relativistic astrometry more than 10 years ago.

Conflicts of Interest

The author declares no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

References

1. Hoffmann-Wellenhof, B.; Moritz, H. Physical Geodesy; Springer: Wien, Austria, New York, NY, USA, 2005.

2. Xie, Y. Relativistic Time Transfer for Inter-satellite Links. Front. Astron. Space Sci.; 2016; 3, 15. [DOI: https://dx.doi.org/10.3389/fspas.2016.00015]

3. Shuai, P. Understanding Pulsars and Space Navigations. Navigation Science and Technology; National Defence Industry Press + Springer Nature Pte Ltd.: Singapore, 2021; Volume 11.

4. Kim, G.; Park, S.; Seong, S.; Lee, J.; Choi, S.; Kim, Y.; Ryu, H.; Lee, S.; Choi, J.; Han, S. Concept of Laser Crosslink Systems Using Nanosatellites in Formation Flying. Acta Astronaut.; 2023; 211, 877. [DOI: https://dx.doi.org/10.1016/j.actaastro.2023.07.014]

5. Zaman, I.; Boyraz, O. Impact of Receiver Architecture on Small Satellite Optical Link in the Presence of Pointing Jitter. Appl. Opt.; 2020; 59, 10177. [DOI: https://dx.doi.org/10.1364/AO.409144]

6. Turyshev, S.; Toth, V.; Sazhin, M. General Relativistic Observables of the GRAIL Mission. Phys. Rev.; 2013; D87, 024020. [DOI: https://dx.doi.org/10.1103/PhysRevD.87.024020]

7. Turyshev, S.; Sazhin, M.; Toth, V. General Relativistic Laser Interferometric Observables of the GRACE-Follow-On Mission. Phys. Rev.; 2014; D89, 105029. [DOI: https://dx.doi.org/10.1103/PhysRevD.89.105029]

8. Turyshev, S.; Yu, N.; Toth, V. General Relativistic Observables for the ACES Experiment. Phys. Rev.; 2016; D93, 045027. [DOI: https://dx.doi.org/10.1103/PhysRevD.93.045027]

9. Atomic Clock Ensemble in Space (ACES). ESA-HSOCOU-020, Rev. 2.0. 2010; Available online: http://wsn.spaceflight.esa.int/docs/Factsheets/20%20ACES%20LR.pdf (accessed on 15 December 2024).

10. Cacciapuoti, L.; Salomon, C. Space Clocks and Fundamental Tests: The ACES Experiment. Eur. Phys. J. Spec. Top.; 2009; 172, pp. 57-68. [DOI: https://dx.doi.org/10.1140/epjst/e2009-01041-7]

11. Shapiro, I. Fourth Test of General Relativity. Phys. Rev. Lett.; 1964; 13, 789. [DOI: https://dx.doi.org/10.1103/PhysRevLett.13.789]

12. Sovers, O.; Fanselow, J.; Jacobs, C. Astrometry and Geodesy with Radio Interferometry: Experiments, Models, Results. Rev. Mod. Phys.; 1998; 70, 1393. [DOI: https://dx.doi.org/10.1103/RevModPhys.70.1393]

13. Gulklett, M. Relativistic Effects in GPS and LEO. Ph.D. Thesis; University of Copenhagen: Copenhagen, Denmark, 8 October 2003; Available online: https://www.yumpu.com/en/document/view/4706552/ (accessed on 30 November 2012).

14. Duchayne, L. Transfert de Temps de Haute Performance: Le Lien Micro-Onde de la Mission ACES. Physique Mathematique. Ph.D. Thesis; Observatoire de Paris: Paris, France, 2008; Available online: https://theses.hal.science/tel-00349882 (accessed on 12 November 2024).

15. Straumann, N. General Relativity; 2nd ed. Springer Science+Business Media: Dordrecht, The Netherlands, 2013.

16. Petit, G.; Wolf, P.; Delva, P. Atomic Time Clocks and Clock Comparisons in Relativistic Spacetimes: A Review. Frontiers in Relativistic Celestial Mechanics; Applications and Experiments De Gruyter and Co.: Berlin, Germany, 2014; Volume 2, 249.

17. Gurfil, P.; Seidelman, P. Celestial Mechanics and Astrodynamics: Theory and Practice; Astrophysics and Space Science Library 436 Springer: Berlin/Heidelberg, Germany, 2016.

18. Kopeikin, S.; Efroimsky, M.; Kaplan, G. Relativistic Celestial Mechanics of the Solar System; Wiley-VCH: New York, NY, USA, 2011.

19. Fock, V. Theory of Space, Time and Gravitation; Pergamonn Press: New York, NY, USA, 1959.

20. Dimitrov, B. Two Null Gravitational Cones in the Theory of GPS-Intersatellite Communications Between Two Moving Satellites. I. Physical and Mathematical Theory of the Space-Time Interval and the Geodesic Distance on Intersecting Null Cones. arXiv; 2017; arXiv: 1712.01101v3

21. Dimitrov, B. Elliptic Integrals for Calculation of the Propagation Time of a Signal, Emitted and Received by Satellites on One Orbit. AIP Conf. Proc.; 2022; 2522, pp. 090004-1-090004-13.

22. Dimitrov, B. Algebraic Geometry and Elliptic Integral Approach for Calculation of the Propagation Time of a Signal Between GPS Satellites with Account of General Relativity Theory. J. Phys. Conf. Ser.; 2023; 2675, 012026. [DOI: https://dx.doi.org/10.1088/1742-6596/2675/1/012026]

23. Pascual-Sanchez, J.-F. Introducing Relativity in Global Navigation Satellite System. Ann. Phys.; 1997; 16, pp. 258-273. [DOI: https://dx.doi.org/10.1002/andp.20075190403]

24. Mishtenko, A.; Fomenko, A. A Course of Differential Geometry and Topology; Factorial Press: Moscow, Russia, 2000.

25. Dimitrov, B. Topological, Differential Geometry Methods and Modified Variational Approach for Calculation of the Propagation Time of a Signal, Emitted by a GPS-Satellite and Depending on the Full Set of 6 Kepler Parameters. J. Phys. Conf. Ser.; 2024; 2910, 012001. [DOI: https://dx.doi.org/10.1088/1742-6596/2910/1/012001]

26. Prasolov, V. Differential Geometry; Springer Nature AG: Cham, Switzerland, 2022.

27. Sivukhin, D. A General Course of Physics, Volume 4 Optics; Fizmatlit Publishing: Moscow, Russia, 2002.

28. Orbital Elements, from Wikipedia, the Free Encyclopedia. Available online: https://en.wikipedia.org/wiki/Orbital_elements (accessed on 2 March 2023).

29. Montenbruck, O.; Gill, E. Satellite Orbits. Models, Methods and Applications; Springer: Berlin/Heidelberg, Germany, 2000.

30. Capderou, M. Handbook of Satellite Orbits: From Kepler to GPS; Springer International Publishing: Cham, Switzerland, 2014.

31. Brauer, D.; Clemence, G. Methods of Celestial Mechanics; Academic Press: New York, NY, USA, London, UK, 1961; pp. 43-60.

32. Xu, G.; Xu, J. Orbits: 2nd Order Singularity-Free Solutions; 2nd ed. Springer: Berlin/Heidelberg, Germany, 2013.

33. Dimitrov, B. New Methods for Analytical Calculation of Elliptic Integrals Applied in Various Physical Problems. AIP Conf. Proc.; 2023; 2953, 070002.

34. Prasolov, V.; Solov’yev, Y. Elliptic Functions and Elliptic Integrals; AMS Translations of Mathematical Monographs 170 AMC: Providence, RI, USA, 1997.

35. Whittaker, E.T.; Watson, G.N. A Course of Modern Analysis: An Introduction; 4th ed. Cambridge University Press: Cambridge, UK, 1927.

36. Nelson, R. Relativistic Time Transfer in the Vicinity of the Earth and in the Solar System. Metrologia; 2011; 48, S171. [DOI: https://dx.doi.org/10.1088/0026-1394/48/4/S07]

37. Ashby, N. Relativity and the Global Positioning System. Phys. Today; 2002; 55, pp. 41-47. [DOI: https://dx.doi.org/10.1063/1.1485583]

38. Fichtenholz, G. Course in Differential and Integral Calculus; 8th ed. Fizmatlit Publishing House: Moscow, Russia, 2003; Volume 2.

39. D’Urso, B.; Etienne, Z.B.; Feldman, M.R.; Kopeikin, S.M.; Trimble, V.; Yi, L. Space-Based Measurements of G. Bull. Am. Astron. Soc.; 2019; 51, 556.

40. Dimitrov, B. New Mathematical Models of GPS Inter-satellite Communications in the Gravitational Field of the Near-Earth Space. AIP Conf. Proc.; 2019; 2075, 040007.

41. Liashko, I.; Boyarchuk, A.; Gai, I.; Kalaida, A. Mathematical Analysis; Publishing House “Vishta Shkola”: Kiev, Ukraine, 1983; Volume 1.

42. Armitage, J.; Eberlein, W. Elliptic Functions; Cambridge University Press: Cambridge, UK, 2006.

43. Mladenov, I.M.; Hadzhilazova, M. The Many Faces of Elastica Forum for Interdisciplinary Mathematics 3; Springer International Publishing AG: Cham, Switzerland, 2017.

44. Legendre, M. Traite des Fonctions Elliptiques es des Integrales Euleriennes. T. 1 Imprimerie de Huzard-Courtier. 1825; Available online: https://gallica.bnf.fr/ark:/12148/bpt6k110147r/f23.item (accessed on 17 October 2007). (In French)

45. Smirnov, V. A Course of Higher Mathematics: Volume 3, Part 2: Complex Variables Special Functions (International Series of Monographs in Pure and Applied Mathematics: Volume 60); Pergamon Press: Oxford, UK, 1964.

46. Timofeev, A. Integration of Functions; State Publishing House for Technical-Theoretical Literature: Moscow, Russia, Leningrad, Russia, 1948.

47. Online Programm. web2.0 Scientific Calculator. Available online: https://web2.0calc.com/ (accessed on 1 September 2024).

48. Murray, C.; Dermott, S. Solar System Dynamics; Cambridge University Press: Cambridge, UK, 2009.

49. Ashby, N. Relativity in GNSS. Springer Handbook of Spacetime; Ashtekar, A.; Petkov, S. Springer: Berlin/Heidelberg, Germany, London, UK, New York, NY, USA, 2014; pp. 509-526.

50. Petit, G.; Wolf, P. Relativistic Theory for Picosecond Time Transfer in the Vicinity of the Earth. Astron. Astrophys.; 1994; 286, 971.

51. Ashby, N. Relativity in the Global Positioning System. Living Rev. Relativ.; 2003; 6, 1. [DOI: https://dx.doi.org/10.12942/lrr-2003-1] [PubMed: https://www.ncbi.nlm.nih.gov/pubmed/28163638]

52. Ashby, N. Relativistic effects in the Global Positioning System. Gravitation and Relativity at the Turn of the Millenium, Proceedings of the 15th International Conference on General Relativity and Gravitation, Pune, India, 16–21 December 1997; Dadhich, N.; Narlikar, J. International University Centre for Astronomy and Astrophysics: Pune, India, 1998.

53. Beard, R.; Senior, R. Clocks. Springer Handbook of Global Navigation Satellite Systems; Teunissen, P.; Montenbruck, O. Springer International Publishing AG: Cham, Switzerland, 2017; pp. 180-190.

54. Kraniotis, G.V. Precise Relativistic Orbits in Kerr and Kerr-(anti) de Sitter Spacetimes. Class. Quantum Gravity; 2004; 21, pp. 4743-4769. [DOI: https://dx.doi.org/10.1088/0264-9381/21/19/016]

55. Husemoller, D. Elliptic Curves. Graduate Texts in Mathematics 111; Springer: New York, NY, USA, Berlin/Heidelberg, Germany, 1987.

56. Lang, S. Elliptic Functions. Graduate Texts in Mathematics 112; 2nd ed. Springer: Berlin/Heidelberg, Germany, New York, NY, USA, 1987.

57. Lavrentiev, M.; Shabat, B. Methods of the Theory of Functions of Complex Variable; 4th ed. Science Publishing House: Moscow, Russia, 1978.

58. Markushevich, A. A Short Course on the Theory of Analytical Functions; Nauka Publishing House: Moscow, Russia, 1966.

59. Knapp, A. Elliptic Curves Mathematical Notes 40; Princeton University Press: Princeton, NJ, USA, 1992; Available online: https://cjhb.site/Files.php/Books/Differential%20Geometry/Knapp.-.Elliptic.Curves%2C.Princeton.1992.223s.pdf (accessed on 1 April 2023).

60. Koblitz, N. Introduction to Elliptic Curves and Modular Forms Graduate Texts in Mathematics 97; Springer: New York, NY, USA, Berlin/Heidelberg, Germany, 1984.

61. Obreshkoff, N. Higher Algebra; State Publishing House “Science and Arts”: Sofia, Bulgaria, 1954.

62. Schur, I. Uber Potenzreihen, die im Innern des EinheitskreisesBeschrankt sind (Fortsetzung). J. Fur Die Reine Und Angewandte Math.; 1918; 148, pp. 122-145. [DOI: https://dx.doi.org/10.1515/crll.1918.148.122]

63. Markushevich, A. Theory of Functions of a Complex Variable; Chelsea: New York, NY, USA, 1977; Volume 1 and 2.

64. Silverman, J. The Arithmetic of Elliptic Curves Graduate Texts in Mathematics 106; 2nd ed. Springer Science+Business Media LLC: Dordrecht, The Netherlands, Berlin/Heidelberg, Germany, 2009.

65. Silverman, J.H.; Tate, J.T. Rational Points on Elliptic Curves Undergraduate Text in Mathematics; 2nd ed. Springer International Publishing: Cham, Switzerland, Berlin/Heidelberg, Germany, New York, NY, USA, 2015.

66. Whittaker, E. A Treatise on the Analytical Dynamics of Particles and Rigid Bodies; Cambridge University Press: Cambridge, UK, 1947.

67. Hancock, H. Elliptic Integrals; John Wiley and Sons: New York, NY, USA, Chapman Hall: London, UK, 1917.

68. Kurosh, A. Higher Algebra; Translation from Russian into English by G. Yankovsky Mir Publishers: Moscow, Russia, 1972; Available online: https://www.math.uni-bielefeld.de/~grigor/kurosh-higher-algebra.pdf (accessed on 3 April 2022).

69. Prasolov, V.V. Polynomials Algorithms and Computation in Mathematics Volume 11; Springer: Berlin/Heidelberg, Germany, 2004.

70. Tabor, M. Chaos and Integrability in Nonlinear Dynamics: An Introduction; John Wiley and Sons Inc.: New York, NY, USA, 1989.