The solution of the Kepler equation is computational

Question

I am trying to solve Kepler Equation as a step towards finding the true anomaly of the orbital body in time . It turns out, however, that the Kepler equation is difficult to solve, and the wikipedia page describes a process using calculus . Well, I don’t know calculus, but I understand that solving an equation involves an infinite number of sets that give a closer approximation to the correct answer.

I don’t see, to look at the math, how to do it computationally, so I was hoping that someone with a better math background could help me. How can I solve this beast computationally?

FWIW, I use F # - and I can calculate the other elements needed for this equation, this is only this part that I came across.

I am also open to methods that approximate the true anomaly of a given time, the distance between periapsis and eccentricity

Answer 1

You can check it out implemented in C # by Carl Johansen

Represents a body in elliptical orbit about a massive central body

Here is a comment from the code

The true anomaly in this context is the angle between the body and the sun. For elliptical orbits, this is a bit complicated. The percentage of the period is still key, but we also need to apply the Kepler equation (based on eccentricity) to ensure that we areas in equal terms. This equation is transcendental (i.e., cannot be solved algebraically), so we either use an approximating equation or solve it using a numerical method. In my implementation, the Newton-Raphson Iteration to get an excellent approximate answer (usually in 2 or 3 iterations).

Answer 2

This article:

Practical Method for Solving the Kepler Equation
http://murison.alpheratz.net/dynamics/twobody/KeplerIterations_summary.pdf

shows how to solve the Kepler equation using an iterative computational method. It should be simple enough to translate it into the language of your choice.

---

You can also find this interesting. This is ocaml, part of which claims to contain a solver for the Kepler equation. Since F # is in the ML language family (like oamaml), this can be a good starting point.

Answer 3

wanted to drop the answer here if this page is found by anyone looking for similar content.

The following was written as an “expression” in Adobe After Effects software, so it is javascriptish, although I have a Python version for another application (movie 4d). The idea is the same: execute the Newton method iteratively until any arbitrary precision is achieved.

Please note that I do not place this code as exemplary or significantly effective in any case, just sending the code that we created as a deadline to perform a specific task (namely, to move the planet around the focus in accordance with Kepler’s laws and do it exactly ) We do not write code for life, and therefore we also do not publish it for criticism. Quick and dirty, on time.

In After Effects, any expression code is executed once - for each frame in the animation. This limits what can be done in the implementation of many algorithms due to the inability to easily cope with global data (other algorithms for Kepler motion use intermediate updated velocity vectors, an approach that we could not use). The result left by the code is the [x, y] position of the object at that moment in time (inside, this is the frame number), and the code is designed to bind the position of the object layer to the element on the chart.

This code developed from material found at http://www.jgiesen.de/kepler/kepler.html , and is suggested here for the next guy.

 pi = Math.PI; function EccAnom(ec,am,dp,_maxiter) { // ec=eccentricity, am=mean anomaly, // dp=number of decimal places pi=Math.PI; i=0; delta=Math.pow(10,-dp); var E, F; // some attempt to optimize prediction if (ec<0.8) { E=am; } else { E= am + Math.sin(am); } F = E - ec*Math.sin(E) - am; while ((Math.abs(F)>delta) && (i<_maxiter)) { E = E - F/(1.0-(ec* Math.cos(E) )); F = E - ec * Math.sin(E) - am; i = i + 1; } return Math.round(E*Math.pow(10,dp))/Math.pow(10,dp); } function TrueAnom(ec,E,dp) { S=Math.sin(E); C=Math.cos(E); fak=Math.sqrt(1.0-ec^2); phi = 2.0 * Math.atan(Math.sqrt((1.0+ec)/(1.0-ec))*Math.tan(E/2.0)); return Math.round(phi*Math.pow(10,dp))/Math.pow(10,dp); } function MeanAnom(time,_period) { curr_frame = timeToFrames(time); if (curr_frame <= _period) { frames_done = curr_frame; if (frames_done < 1) frames_done = 1; } else { frames_done = curr_frame % _period; } _fractime = (frames_done * 1.0 ) / _period; mean_temp = (2.0*Math.PI) * (-1.0 * _fractime); return mean_temp; } //============================== // a=semimajor axis, ec=eccentricity, E=eccentric anomaly // delta = delta digits to exit, period = per., in frames //---------------------------------------------------------- _eccen = 0.9; _delta = 14; _maxiter = 1000; _period = 300; _semi_a = 70.0; _semi_b = _semi_a * Math.sqrt(1.0-_eccen^2); _meananom = MeanAnom(time,_period); _eccentricanomaly = EccAnom(_eccen,_meananom,_delta,_maxiter); _trueanomaly = TrueAnom(_eccen,_eccentricanomaly,_delta); r = _semi_a * (1.0 - _eccen^2) / (1.0 + (_eccen*Math.cos(_trueanomaly))); x = r * Math.cos(_trueanomaly); y = r * Math.sin(_trueanomaly); _foc=_semi_a*_eccen; [1460+x+_foc,540+y];