Kepler's Equation Research Articles

Machine learning meets Kepler: inverting Kepler’s equation for All vs All conjunction analysis

The number of satellites in orbit around Earth is increasing rapidly, with the risk of collision rising accordingly. Trends of the global population of satellites need to be analyzed to test the viability and impact of proposed rules and laws affecting the satellite population and collision avoidance strategies. This requires large scale simulations of satellites that are propagated on long timescales to compute the large amounts of actionable close encounters (called conjunctions), which could lead to collisions. Rigorously checking for conjunctions by computing future states of orbits is computationally expensive due to the large amount of objects involved and conjunction filters are thus used to remove non-conjuncting orbit pairs from the list of possible conjunctions. In this work, we explore the possibility of machine learning (ML) based conjunction filters using several algorithms such as eXtreme Gradient Boosting, TabNet and (physics-informed) neural networks and deep operator networks. To show the viability and the potential of ML based filters, these algorithms are trained to predict the future state of orbits. For the physics-informed approaches, multiple partial differential equations are set up using the Kepler equation as a basis. The empirical results demonstrate that physics-informed deep operator networks are capable of predicting the future state of orbits using these equations (RMSE: 0.136) and outperform eXtreme Gradient Boosting (RMSE: 0.568) and TabNet (RMSE: 0.459). We also propose a filter based on the trained deep operator network which is shown to outperforms the filter capability of the commonly used perigee-apogee test and the orbit path filter on a synthetic dataset, while being on average 3.2 times faster to compute than a rigorous conjunction check.

Finding the Exact Solution of Kepler’s Equation for an Elliptical Satellite Orbit Using the First Kind Bessel Function

In this study, the first kind Bessel function was used to solve Kepler equation for an elliptical orbiting satellite. It is a classical method that gives a direct solution for calculation of the eccentric anomaly. It was solved for one period from (M=0-360)° with an eccentricity of (e=0-1) and the number of terms from (N=1-10). Also, the error in the representation of the first kind Bessel function was calculated. The results indicated that for eccentricity of (0.1-0.4) and (N = 1-10), the values of eccentric anomaly gave a good result as compared with the exact solution. Besides, the obtained eccentric anomaly values were unaffected by increasing the number of terms (N = 6-10) for eccentricities (0.8 and 0.9). The Bessel function's solution appeared to be close to the exact solution for eccentricity of 1 and more than 10 number of terms. Finally, the representation of the first kind Bessel function J1(x) was closer to the exact representation only for eccentricity 0.5 and (N=1-10).

On the Bessel Solution of Kepler’s Equation

Since its introduction in 1650, Kepler’s equation has never ceased to fascinate mathematicians, scientists, and engineers. Over the course of five centuries, a large number of different solution strategies have been devised and implemented. Among them, the one originally proposed by J. L. Lagrange and later by F. W. Bessel still continue to be a source of mathematical treasures. Here, the Bessel solution of the elliptic Kepler equation is explored from a new perspective offered by the theory of the Stieltjes series. In particular, it has been proven that a complex Kapteyn series obtained directly by the Bessel expansion is a Stieltjes series. This mathematical result, to the best of our knowledge, is a new integral representation of the KE solution. Some considerations on possible extensions of our results to more general classes of the Kapteyn series are also presented.

A new method for solving the hyperbolic Kepler equation

This paper proposes a highly efficient method for solving the hyperbolic Kepler equation (HKE). The hyperbolic eccentric anomaly interval is divided into two parts: a finite interval and an infinite interval. For the finite interval, a piecewise Padé approximation is first used to establish an initial approximate solution of the HKE. For the infinite interval, an analytical initial approximate solution of the HKE is constructed. These initial approximations are highly accurate and can be further improved to higher accuracy with only one step of Schröder iteration. The proposed method only requires the evaluation of no more than three transcendental functions.

A functionally-fitted block hybrid Falkner method for Kepler equations and related problems

A functionally-fitted block hybrid Falkner method for Kepler equations and related problems

Orbit Predictions Using KS Elements With Earth’s Flattening

Predicting the orbit of a satellite is a fundamental requirement in many areas of aerospace, such as mission planning, satellite geodesy, re-entry prediction, maneuver planning, collision avoidance, formation flying, etc. For near-Earth orbits, the forces due to the non spherical nature of the Earth and atmospheric drag plays an important role and are mainly responsible to bring the satellite back to the Earth. Thus inclusion of the effect of these perturbing forces becomes important for precise orbit computation of near-Earth orbits. Many transformations have emerged in the literature to stabilize the equations of motion either to reduce the accumulation of local numerical errors or allowing the use of large integration step sizes, or both in the transformed space. One such transformation is known as KS transformation by Kustaanheimo and Stiefel [1], who regularized the nonlinear Kepler equation of motion and reduced it into linear differential equations of a harmonic oscillator of constant frequency. In this paper a detailed analysis is carried out for orbit prediction using KS differential equations by including the non spherical gravitational potential of the Earth as the perturbing force. Earths zonal and tesseral harmonics, which produce bands of constant deviation from the spherical field along lines of latitude and longitude, are included to precisely model the Earth’s gravitational potential. Higher order Earth’s flattening terms are included by utilizing the recurrence relations of associated Legendre polynomial and its derivatives. To know the effectiveness of the theory, the results are compared with the real IRS-1A satellite data and the results of other existing solutions for a long duration of nearly a month time. The comparison shows that the KS method provides one of the most and accurate techniques for orbit prediction available at present.

Refining the OJ 287 2022 impact flare arrival epoch

ABSTRACT The bright blazar OJ 287 routinely parades high brightness bremsstrahlung flares, which are explained as being a result of a secondary supermassive black hole (SMBH) impacting the accretion disc of a more massive primary SMBH in a binary system. The accretion disc is not rigid but rather bends in a calculable way due to the tidal influence of the secondary. Next, we refer to this phenomenon as a variable disc level. We begin by showing that these flares occur at times predicted by a simple analytical formula, based on general relativity inspired modified Kepler equation, which explains impact flares since 1888. The 2022 impact flare, namely flare number 26, is rather peculiar as it breaks the typical pattern of two impact flares per 12-yr cycle. This is the third bremsstrahlung flare of the current cycle that follows the already observed 2015 and 2019 impact flares from OJ 287. It turns out that the arrival epoch of flare number 26 is sensitive to the level of primary SMBH’s accretion disc relative to its mean level in our model. We incorporate these tidally induced changes in the level of the accretion disc to infer that the thermal flare should have occurred during 2022 July–August, when it was not possible to observe it from the Earth. Thereafter, we explore possible observational evidence for certain pre-flare activity by employing spectral and polarimetric data from our campaigns in 2004/05 and 2021/22. We point out theoretical and observational implications of two observed mini-flares during 2022 January–February.

Calculation Mars – Earth distance and Mars orbital elements with Julian date

In this paper, the Mars orbital elements were calculated. These orbital elements—the major axis, the inclination (i), the longitude of the ascending node (W), the argument of the perigee (w), and the eccentricity (e)—are essential to knowing the size and shape of Mars' orbit. The quick basic program was used to calculate the orbital elements and distance of Mars from the Earth from 25/5/1950 over 10000 days. These were calculated using the empirical formula of Meeus, which depended on the Julian date, which slightly changed for 10000 days; Kepler's equation was solved to find Mars' position and its distance from the Sun. The ecliptic and equatorial coordinates of Mars were calculated. The distance between Mars and the center of the Earth, in astronomical units (A.U.), was calculated. RM-E(min) was found to be between 0.4763 and 0.5108, and RM-E(max) was found to be between 2.548 and 2.6259. Furthermore, the findings revealed that the Mars orbital elements have changed over time.

Analytical images of Kepler's equation solutions and their applications

The simple fast-converging analytical calculations algorithms for eccentric anomaly are proposed for an arbitrary eccentricity 0<e≤1. The kinematic characteristics of Halley's comet are calculated as the function of time. Mass of Galaxy + NGC 224 system using the model of elliptical relative motion is also estimated.

A new solution approach via analytical approximation of the elliptic Kepler equation

A new analytical approach for constructing approximate solutions to the elliptic Kepler equation is proposed. We first establish a high-accuracy initial approximation using the piecewise Padé approximation, subsequently we apply the Schröder method to further improve the accuracy of the initial approximation. In general, one Schröder iteration is sufficient to obtain a highly accurate approximate solution. This is a direct method that requires only solving a cubic equation and evaluating two trigonometric functions. The approximate, analytical solutions are compared with solutions by other numerical procedures to prove the accuracy and effectiveness of the proposed approach.

Solar Position Algorithm Based on the Kepler Equation

When calculating the position of the sun, earth’s motion can be assumed to be an ellipse if the accuracy of calculation is required to be 0.01 degrees. Then, Kepler’s equation can be applied from the mean anomaly of the sun at a specific time to calculate the true anatomy of the sun at the time, and the sun’s position can be calculated. The average absolute error of calculating the sun’s altitude and azimuth is only 0.04 and 0.06 mrad, respectively, which can meet the requirements of a concentrated solar tracking system. This method only needs to correct the length of the regression year and the near point year for every 100 years, so it can be used for a long time.

Two fast and accurate routines for solving the elliptic Kepler equation for all values of the eccentricity and mean anomaly

Context. The repetitive solution of Kepler’s equation (KE) is the slowest step for several highly demanding computational tasks in astrophysics. Moreover, a recent work demonstrated that the current solvers face an accuracy limit that becomes particularly stringent for high eccentricity orbits. Aims. Here we describe two routines, ENRKE and ENP5KE, for solving KE with both high speed and optimal accuracy, circumventing the abovementioned limit by avoiding the use of derivatives for the critical values of the eccentricity e and mean anomaly M, namely e > 0.99 and M close to the periapsis within 0.0045 rad. Methods. The ENRKE routine enhances the Newton-Raphson algorithm with a conditional switch to the bisection algorithm in the critical region, an efficient stopping condition, a rational first guess, and one fourth-order iteration. The ENP5KE routine uses a class of infinite series solutions of KE to build an optimized piecewise quintic polynomial, also enhanced with a conditional switch for close bracketing and bisection in the critical region. High-performance Cython routines are provided that implement these methods, with the option of utilizing parallel execution. Results. These routines outperform other solvers for KE both in accuracy and speed. They solve KE for every e ∈ [0, 1 − ϵ], where ϵ is the machine epsilon, and for every M, at the best accuracy that can be obtained in a given M interval. In particular, since the ENP5KE routine does not involve any transcendental function evaluation in its generation phase, besides a minimum amount in the critical region, it outperforms any other KE solver, including the ENRKE, when the solution E(M) is required for a large number N of values of M. Conclusions. The ENRKE routine can be recommended as a general purpose solver for KE, and the ENP5KE can be the best choice in the large N regime.

New methods of isochrone mechanics

Isochrone potentials are spherically symmetric potentials within which a particle orbits with a radial period that is independent of its angular momentum. Whereas all previous results on isochrone mechanics have been established using classical analysis and geometry, in this article, we revisit the isochrone problem of motion using tools from Hamiltonian dynamical systems. In particular, we (1) solve the problem of motion using a well-adapted set of angle-action coordinates and generalize the notion of eccentric anomaly to all isochrone orbits, and (2) we construct the Birkhoff normal form for a particle orbiting a generic radial potential and examine its Birkhoff invariants to prove that the class of isochrone potentials is in correspondence with parabolas in the plane. Along the way, several fundamental results of celestial mechanics, such as the Bertrand theorem or the Kepler equation and laws, are obtained as special cases of more general properties characterizing isochrone mechanics.

Comment on ‘An efficient code to solve the Kepler equation: elliptic case’

ABSTRACTIn a recent MNRAS article, Raposo-Pulido and Pelaez (RPP) designed a scheme for obtaining very close seeds for solving the elliptic Kepler equation with the classical and modified Newton–Raphson methods. This implied an important reduction in the number of iterations needed to reach a given accuracy. However, RPP also made strong claims about the errors of their method that are incorrect. In particular, they claim that their accuracy can always reach the level of ∼5ε, where ε is the machine epsilon (e.g. ε = 2.2 × 10−16 in double precision), and that this result is attained for all values of the eccentricity e < 1 and the mean anomaly M ∈ [0, π], including for e and M that are arbitrarily close to 1 and 0, respectively. However, we demonstrate both numerically and analytically that any implementation of the classical or modified Newton–Raphson methods for Kepler’s equation, including those described by RPP, has a limiting accuracy of the order of ${\sim}\varepsilon /\sqrt{2(1-e)}$. Therefore the errors of these implementations diverge in the limit e → 1, and differ dramatically from the incorrect results given by RPP. Despite these shortcomings, the RPP method can provide a very efficient option for reaching such limiting accuracy. We also provide a limit that is valid for the accuracy of any algorithm for solving Kepler equation, including schemes like bisection that do not use derivatives. Moreover, similar results are also demonstrated for the hyperbolic Kepler equation. The methods described in this work can provide guidelines for designing more accurate solutions of the elliptic and hyperbolic Kepler equations.

Bivariate Infinite Series Solution of Kepler’s Equations

A class of bivariate infinite series solutions of the elliptic and hyperbolic Kepler equations is described, adding to the handful of 1-D series that have been found throughout the centuries. This result is based on an iterative procedure for the analytical computation of all the higher-order partial derivatives of the eccentric anomaly with respect to the eccentricity e and mean anomaly M in a given base point (ec,Mc) of the (e,M) plane. Explicit examples of such bivariate infinite series are provided, corresponding to different choices of (ec,Mc), and their convergence is studied numerically. In particular, the polynomials that are obtained by truncating the infinite series up to the fifth degree reach high levels of accuracy in significantly large regions of the parameter space (e,M). Besides their theoretical interest, these series can be used for designing 2-D spline numerical algorithms for efficiently solving Kepler’s equations for all values of the eccentricity and mean anomaly.

An Efficient Iterative Method for Solving the Elliptical Kepler’s Equation

In this paper, a numerical technique, based on Banach’s fixed point theorem, is proposed to obtain an approximate solution of the elliptical Kepler equation which is often used to describe the motion of a planet orbiting its sun. The rapid convergence of the method results in qualitatively accurate solutions in relatively few iterations. In fact, the estimate of approximation error is given as a function of the number of steps of the method. Our results demonstrate that the presented method is very effective and reliable, the solution exists and is unique, provided that the nonlinear term is a contraction, therefore can be considered a very useful and valuable method.

Francesco Carlini: Kepler's equation and the asymptotic solution to singular differential equations

Carlini's career was mainly dedicated to astronomy, but he was also a particularly skilled mathematician. In this article we collect and analyse his mathematical contributions in detail. In particular, in his important Memoir of the year 1817 devoted to Kepler's equation he introduced an innovative idea to solve ordinary differential equations with singular perturbations by means of asymptotic expansions. In the same Memoir also appeared, five years before Laplace's contributions, what is usually called the Laplace limit constant. Furthermore, Carlini published other mathematical Memoirs anticipating, 70 years in advance, the importance of complex branches of the Lambert's special function.

Integration of Deprit's radial intermediary

In Artificial Satellite Theory, for the so called main problem (the two-body problem is only perturbed by the equatorial planet's bulge) Deprit's intermediary plays an important role because of its simplicity and integrability while maintaining a similar structure to the original Hamiltonian. To integrate it, we prove that there is no need to introduce a generalizad Kepler's equation, but the classical Kepler's equation, which has been extensively studied for more than three centuries.

On the Minimality of Keplerian Arcs with Fixed Negative Energy

We revisit a classical result by Jacobi (J Reine Angew Math 17:68–82, 1837) on the local minimality, as critical points of the corresponding energy functional, of fixed-energy solutions of the Kepler equation joining two distinct points with the same distance from the origin. Our proof relies on the Morse index theorem, together with a characterization of the conjugate points as points of geodesic bifurcation.

Re-Evaluation Solution Methods for Kepler's Equation of an Elliptical Orbit

An evaluation was achieved by designing a matlab program to solve Kepler’s equation of an elliptical orbit for methods (Newton-Raphson, Danby, Halley and Mikkola). This involves calculating the Eccentric anomaly (E) from mean anomaly (M=0°-360°) for each step and for different values of eccentricities (e=0.1, 0.3, 0.5, 0.7 and 0.9). The results of E were demonstrated that Newton’s- Raphson Danby’s, Halley’s can be used for e between (0-1). Mikkola’s method can be used for e between (0-0.6).The term that added to Danby’s method to obtain the solution of Kepler’s equation is not influence too much on the value of E. The most appropriate initial Gauss value was also determined to be (En=M), this initial value gave a good result for (E) for these methods regardless the value of e to increasing the accuracy of E. After that the orbital elements converting into state vectors within one orbital period within time 50 second, the results demonstrated that all these four methods can be used in semi-circular orbit, but in case of elliptical orbit Danby’s and Halley’s method use only for e ≤ 0.7, Mikkola’s method for e ≤ 0.01 while Newton-Raphson uses for e < 1, which considers more applicable than others to use in semi-circular and elliptical orbit. The results gave a good agreement as compared with the state vectors of Cartosat-2B satellite that available on Two Line Element (TLE).