Solving Kepler's Equation with Newton-Raphson Method

Question

Note (2022/03/07): This question is solved. Unfortunately, I'm not able to accept the correct answer by Lutz Lehmann, because I screwed up my registration and the account which posted this question is different from the one I'm using now, even though they have the same name.

I'm trying to find the eccentric anomaly of an orbit through the mean anomaly. I need this to calculate the (2D-)position of an orbiter (called Newton orbiter in the following) at different times. I'm using C# and the Unity game engine.

To control the algorithms correctness, I used Unity physics and a simple gravity simulation to simulate another orbiter (called Unity orbiter in the following) with identical properties. Both orbiters completed their orbits in the same time (which was equal to the calculated orbital period and therefore most likely correct). However, they had different velocities along the path. I printed the specific orbital energy of both orbiters each frame. Only the orbital energy of the Unity orbiter remained constant along the orbit.

Since the path and orbital period of my custom orbital calculation are correct and only the velocity along the path is wrong, I suspect that either the mean anomaly or eccentric anomaly calculation is wrong.

The Unity orbiter reaches its top speed correctly at its periapsis (0.7|-3). The Newton orbiter reaches its top speed at (2|6). The gravity well is at (0|0), both bodies start at (4|0) with initial velocity (1|3).

Potential Problems:

• Unity uses a left-handed coordinate system, this might screw up some formulas
• I do not clamp the time and angles to [0, 2*pi], although I believe it should be fine, because sin and cos are periodic

Used Formulas:

Specific Orbital Energy:

$$ \epsilon = \frac {|v|^2} {2} - \frac {\mu} {|r|} $$

Eccentricity Vector:

$$ e = (\frac{|v|^2}{\mu} - \frac{1}{|r|}) * r - \frac{r * v}{\mu} * v $$

Orbital Period:

$$ T = 2 * \pi * \sqrt{\frac{a^3}{\mu}} $$

Mean Anomaly:

$$ M = \frac{t}{T} * 2 * \pi - \phi $$

Newton-Raphson for Eccentric Anomaly:

$$ E_{n+1} = E_n - \frac{f(E_n)}{f'(E_n)} = E_n - \frac{E_n - |e| * sin(E_n) -M}{1 - |e| * cos(E_n)} $$

v - orbiter velocity

$ \mu $ - gravitational parameter (gravity well mass * gravitational constant)

r - orbiter position

a - semimajor axis

t - time since program start in seconds

$ \phi $ - angle between periapsis and x-axis

Code for the Newton orbiter:

public class GravityAnalyzer : MonoBehaviour
{
    public const double GRAVITY_CONSTANT = 0.000000000066743;

    [SerializeField] private double gravityWellMass = 500.0;
    [SerializeField] private Vector2 startVelocity = Vector2.up;
    [SerializeField] private int maxIterations = 10;
    [SerializeField] private double minPrecision = 0.0001;
    private new Rigidbody2D rigidbody = null;
    private double gravitationalParameter = 0.0;
    private double specificOrbitalEnergy = 0.0;
    private double eccentricityMagnitude = 0.0;
    private double semiMajorAxis = 0.0;
    private double semiMinorAxis = 0.0;
    private double phi = 0.0;
    private Vector2 orbitCenter = Vector2.zero;
    private double orbitalPeriod = 0.0;

    Vector2 lastPos = Vector2.zero;

    private void Start()
    {
        rigidbody = gameObject.GetComponent<Rigidbody2D>();

        // https://en.wikipedia.org/wiki/Orbital_elements
        gravitationalParameter = gravityWellMass * GRAVITY_CONSTANT * 1000000000.0; // Dont mind the magic Number, will be cleaned up later
        specificOrbitalEnergy = (startVelocity.sqrMagnitude * 0.5) - (gravitationalParameter / rigidbody.position.magnitude);
        Vector2 eccentricity = ((((startVelocity.sqrMagnitude / (float) gravitationalParameter) - (1.0f / rigidbody.position.magnitude)) * rigidbody.position)
            - ((Vector2.Dot(rigidbody.position, startVelocity) / (float) gravitationalParameter) * startVelocity));
        eccentricityMagnitude = Math.Sqrt((double) eccentricity.x * (double) eccentricity.x + (double) eccentricity.y * (double) eccentricity.y);
        // https://en.wikipedia.org/wiki/Semi-major_and_semi-minor_axes#Energy;_calculation_of_semi-major_axis_from_state_vectors
        semiMajorAxis = -gravitationalParameter / (2.0 * specificOrbitalEnergy);
        semiMinorAxis = Math.Sqrt(1.0 - eccentricityMagnitude * eccentricityMagnitude) * semiMajorAxis;
        // https://en.wikipedia.org/wiki/Orbital_period
        orbitalPeriod = 2.0 * Math.PI * Math.Sqrt((semiMajorAxis * semiMajorAxis * semiMajorAxis) / gravitationalParameter);

        // phi is the Angle by which the Orbit is rotated around the Origin of the Coordinate System
        // Vector2.SignedAngle() returns the Result in Degrees, must convert to Radians!
        phi = Vector2.SignedAngle(eccentricity, Vector2.right) * (Math.PI / 180.0);

        // Calculate the initial Eccentric Anomaly to be able to solve for the Center Point in the next Step
        double eccentricAnomaly = CalculateEccentricAnomaly();

        // Calculate Center Point of the Ellipse
        orbitCenter = -(new Vector2((float) (semiMajorAxis * Math.Cos(phi) * Math.Cos(eccentricAnomaly) - semiMinorAxis * Math.Sin(phi) * Math.Sin(eccentricAnomaly)),
            (float) (semiMajorAxis * Math.Sin(phi) * Math.Cos(eccentricAnomaly) + semiMinorAxis * Math.Cos(phi) * Math.Sin(eccentricAnomaly)))
            - rigidbody.position);
    }

    private void FixedUpdate()
    {
        // Calculate Eccentric Anomaly from current Time
        double eccentricAnomaly = CalculateEccentricAnomaly();

        // Plug in all Parameters in Parameter-Form of Ellipse-Equation
        // https://de.wikipedia.org/wiki/Ellipse#Ellipsengleichung_(Parameterform)
        rigidbody.position = orbitCenter
            + new Vector2((float) (semiMajorAxis * Math.Cos(phi) * Math.Cos(eccentricAnomaly) - semiMinorAxis * Math.Sin(phi) * Math.Sin(eccentricAnomaly)),
            (float) (semiMajorAxis * Math.Sin(phi) * Math.Cos(eccentricAnomaly) + semiMinorAxis * Math.Cos(phi) * Math.Sin(eccentricAnomaly)));
    }

    private double CalculateEccentricAnomaly()
    {
        // Calculate Mean Anomaly
        double meanAnomaly = (Time.time / orbitalPeriod) * 2.0 * Math.PI - phi;

        // Solve Kepler Equation and convert Mean Anomaly to Eccentric Anomaly
        // https://en.wikipedia.org/wiki/Kepler%27s_equation#Numerical_approximation_of_inverse_problem
        double eccentricAnomaly = meanAnomaly;
        // Skip Calculation if Orbit is circular
        if(eccentricityMagnitude > 0.01)
        {
            // Use different initial Guess for very elliptic Orbits
            if(eccentricityMagnitude > 0.8)
            {
                eccentricAnomaly = Math.PI;
            }

            // Newton-Raphson
            int i = 0;
            double lastEccentricAnomaly = 0.0;
            do
            {
                lastEccentricAnomaly = eccentricAnomaly;
                eccentricAnomaly = eccentricAnomaly
                    - ((eccentricAnomaly - eccentricityMagnitude * Math.Sin(eccentricAnomaly) - meanAnomaly)
                    / (1.0 - eccentricityMagnitude * Math.Cos(eccentricAnomaly)));
                Debug.Log(i + ": Last: " + lastEccentricAnomaly + " Current: " + eccentricAnomaly);
            }
            while(i++ < maxIterations && !((eccentricAnomaly - lastEccentricAnomaly) < minPrecision && (eccentricAnomaly - lastEccentricAnomaly) > -minPrecision));
        }

        return eccentricAnomaly;
    }
}

Approximate Orbit Sequence:

• Green is the Newton orbiter, red is the Unity orbiter
• Initial Velocity is (1|3) for both

$$ t = 0 $$

$$ t ~= 0.5 * T $$

$$ t ~= 0.9 * T $$

Answer 1

In summary: it is wrong to use the phase angle of the periapsis also as the offset for the mean anomaly. Both transformations, mean to eccentric anomaly, and eccentric anomaly to true anomaly (and to the angle in the original coordinate system) are non-linear.

Find constants for the Kepler laws

I'll first run my Python "standard program" to compute the orbit geometry. I use the variables

• $a,b,c$ for the long and short semi-axis and the distance from the center to the focal point of the ellipse, $a^2=b^2+c^2$,
• $e$ the eccentricity, $c=e·a$ thus $b=\sqrt{1-e^2}a$
• $\phi$ is the angle in the "world" coordinate system, $\psi=\phi_{\rm peri}$ the angle of the periapsis in world coordinates, $\theta=\phi-\psi$ the "true anomaly", the angle in the rotated coordinate system parallel to the long and short axis, still having the focal point as origin.

import math as m
def polar(x,y): return m.hypot(x,y), m.atan2(y,x)

# input
GM = 0.000000000066743*500*1000000000.0;

x0,y0,vx0,vy0 = 4,0,1,3

# transform initial conditions to polar coordinates
r0, phi0 = polar(x0,y0)
dotr0 = (x0*vx0+y0*vy0)/r0
L = y0*vx0-x0*vy0    # r^2*dotphi = L constant, L^2 = GM*R
dotphi0 = L/r0**2

# find parameters of the ellipse in polar coordinates
ex = -vy0*L/GM-x0/r0; ey = +vx0*L/GM-y0/r0
e, psi = polar(ex,ey)

R = L**2/(GM)

# compute geometry of the ellipse in Cartesian coordinates
a=R/(1-e**2); b=R/(1-e**2)**0.5; c=a*e

# period of the orbit
T = 2*m.pi/(GM)**0.5*a**1.5

print(f"a={a:.6f}, b={b:.6f}, c={c:.6f},") 
print(f"eccentricity e={e:.6f}, periapsis angle psi={psi/m.pi:.6f}*pi")
print(f"periapsis = {a-c:.6f}, apoapsis = {a+c:.6f}")
print(f"orbital period T={T:.6f}")

with the result

a=4.991437, b=4.640942, c=1.837416,
eccentricity e=0.368114, periapsis angle psi=-0.431361*pi
periapsis = 3.1540, apoapsis = 6.82885285
orbital period T=12.129151

As far as I can see, your initialization should produce the same values.

True to eccentric anomaly -- angle to focus vs. angle to center

Now the relation between true anomaly and eccentric anomaly is given by the tangent relation $$ \tan(E/2)=\sqrt{\frac{1-e}{1+e}}\tan(\theta/2) $$ This gives the eccentric anomaly at time $t=0$ as $$ E_0=2\arctan\left(\sqrt{\frac{1-e}{1+e}}\tan(\theta_0/2)\right) $$ where $\theta_0=\phi_0-\phi_{\rm peri}$. In most cases both of the angles are non-trivial.

Eccentric to mean anomaly - angle to center vs. time

From this the mean anomaly can be reconstructed as $$ M_0=2\pi\frac{0-t_{\rm peri}}{T}=E_0-e\sin(E_0). $$ Any other mean anomaly can now be computed as $$ M=2\pi\frac{t-t_{\rm peri}}{T}=M_0+2\pi\frac{t}{T}. $$ Or use $M_{\rm peri}=-M_0$ to get a difference to a base point. Note that you used $\phi_{\rm peri}$ in the place of $M_{\rm peri}$. While both angles fill the full circle with a monotonic relation, and are simultaneously zero, they are not equal, especially for medium to large eccentricity $e$.

Then apply the Newton method to get the eccentric anomaly and from that the position data corresponding to time $t$.

Exact solution vs. numerical solution

Comparing that to the numerical solution direct from the initial conditions can be done by plotting them over each other with with markers and different colors and line widths.

E0 = 2 * m.atan(((1-e)/(1+e))**0.5 * m.tan((phi0-psi)/2))
M0 = E0 - e*m.sin(E0)

# solve M = (1-e)*E + e*(E-sin(E)) for E
def ecccA(t):
    M = M0+2*m.pi*t/T
    E=M
    for _ in range(10):
        dE = -(E-e*m.sin(E)-M)/(1-e*m.cos(E))
        E  = E+dE
        if abs(dE)<1e-7: break;
    return E

print(E0, ecccA(0))

# E to rotated Cartesian frame
def pos(t):
    E = ecccA(t)
    xo = a*m.cos(E)-c
    yo = b*m.sin(E)
    return xo*m.cos(psi)-yo*m.sin(psi), +xo*m.sin(psi)+yo*m.cos(psi)

print([x0,y0],pos(0))

# numerical integration
## gravity equations
def eqn(u,t):
    x,y,vx,vy = u
    r3 = (x*x+y*y)**1.5
    return vx,vy,-GM*x/r3,-GM*y/r3

Refine=10
to = np.linspace(0,T,12*Refine+1)

u = odeint(eqn,[x0,y0,vx0,vy0],to).T

kepler = np.asarray([pos(t) for t in to]).T

plt.figure(figsize = (5,5))
plt.gca().set_aspect('equal')
plt.plot(x0,y0,'go',ms=10)
plt.plot(u[0],u[1],'b', lw=5)
plt.plot(u[0,::Refine],u[1,::Refine],'bo',ms=8)
plt.plot(kepler[0,0],kepler[1,0],'yo',ms=4)
plt.plot(kepler[0],kepler[1],'y',lw=3)
plt.plot(kepler[0,::Refine],kepler[1,::Refine],'y+',ms=6)

plt.grid(); plt.show()

The result now indeed fits perfectly