I am trying to model a two dimensional orbit for a two body Kepler problem but have gotten stuck when introducing the time variable.

For a satellite with known semi major axis (a), eccentricity (e), and true anomaly (theta), I have:

r = a*(1-e**2)/(1+e*cos(theta))

How can I calculate theta as a function of time using the Wiki suggested Universal variable formulation method? I have no idea how to implement (am using Python but any algorithm advice much appreciated!)

Alternatively, how do I calculate r as a function of time?

Note: All other orbital elements and masses are available. Also I am trying to come up with a general solution for elliptical, hyperbolic and parabolic orbits.


  • \$\begingroup\$ I see that you used the tags game-design and game-mechanics. Is this closed form really necessary for a game-level simulation, and how does it interact with the gameplay? \$\endgroup\$ Dec 3 '15 at 8:27
  • \$\begingroup\$ As for your problem, you might be interested in this book which targets just an audience like you (and me): amazon.com/Astronomical-Algorithms-Jean-Meeus/dp/0943396352 \$\endgroup\$ Dec 3 '15 at 10:02
  • \$\begingroup\$ Lars- I'm opting for closed form (if I have understood what that means correctly) to give planets fixed, stable orbits over a long period of time. Smaller object (spaceships etc) will have variable motion defined by newtons gravitational laws, which should be easy enough to calculate. \$\endgroup\$
    – imnegan
    Dec 3 '15 at 11:02
  • \$\begingroup\$ Planetmaker- thanks for the recommendation, I'll check it out. \$\endgroup\$
    – imnegan
    Dec 3 '15 at 11:03

I calculated true anomaly as function of time, for planetary motion , in c# , in this way:

  1. Compute mean anomaly (time: current time , G: newton grav.connstant, M: planet mass or the sum of the two orbiting objects , a: semi major axis)

            //M = nt
                double n = Math.Sqrt((G * (M)) / (a * a * a));
                double Mt = n * time;
  2. Compute the eccentric anomaly E by solving Kepler's equation:

         //For orbits with ε > 0.8, an initial value of E0 = π should be used.
            if (eccentr>0.8)
                E = NumApprox(150, Math.PI,Mt, 10E-15);
                E = NumApprox(150, Mt, Mt, 10E-15);
  3. true anomaly (angle)

    true_anom = 2.0 * Math.Atan2(Math.Sqrt(1.0 + eccentr) * Math.Sin(E / 2.0), Math.Sqrt(1.0 - eccentr) * Math.Cos(E / 2.0));

  4. distance from planet

    d = a * ((1.0 - eccentr * eccentr) / (1.0 + eccentr * Math.Cos(true_anom)));

Finaly Numerical approximation of inverse problem:

private double NumApprox(int intr, double prev,double Mt, double  err)
        double ret = prev;
        double retprev = prev;
        for (int i=0 ; i<intr; i++){
            retprev = ret;
            ret = ret - (ret - eccentr * Math.Sin(ret) - Mt) / (1.0 - eccentr * Math.Cos(ret));
            if ( Math.Abs(ret - retprev) < err)
        return ret; 

enter image description here

EDIT: calculate Position and Velocity

What we done is half of the work in getting Cartesain Orbit Elements from kepler Orbit Elements where :

    //some kepler Orbit Elements: 
    public double d;
    public double true_anom;
    public double eccentr;
    public double a;
    public double E;
    public double w=0; //small omega ω :  Argument of periapsis (in rad)

    //Cartesain Orbit Elements: 
    //Position Vector
    private double x;
    private double y;
    private double z; //2d :not used
    //Velocity Vector
    private double vx;
    private double vy;
    private double vz; //2d :not used

as final step we can calculate position and velocity vectors as :

x = d * Math.Cos(true_anom);
y = d * Math.Sin(true_anom);
//apply ω
double xx = x * Math.Cos(w) - y * Math.Sin(w);
double yy = x * Math.Sin(w) + y * Math.Cos(w);
x = xx;
y = yy;

double v = Math.Sqrt(G * M * a) / d;
vx = -v * Math.Sin(E);
vy = -v * Math.Sqrt(1.0-eccentr*eccentr) * Math.Cos(E);

EDIT: references enter image description here

Orbital elements

Kepler's equation

True anomaly

Eccentricity vector

Argument of periapsis

Newton's method for numerical approximation

The Kerbal Space Program (KSP) Physics Documantation (pdf)

  • \$\begingroup\$ This is great, thanks! Also the fact that it's in 2D just happens to align with what I'm toying with. Out of curiosity, without using inclination and the z axis is it possible to represent objects orbiting in a clockwise direction? \$\endgroup\$
    – imnegan
    Apr 17 '19 at 12:17

Retro orbit function, also can be adapted to static-node elements (ap,ma) ... this is the math. result error +- 0.08 degree from VSOP.

function skyUpdate3(id) { 
/// Using timed Elements - delivered in radians.
var e = id.e, n, E, r, v, o, x = [0,0,0];
var ap = id.lp - id.o;
var M = toNormalRange(id.ml - id.lp);
  E = SolveKepler(M,e);
  r = IDS(au2km(id.a)) * (1.0 - (e * Math.cos(E)));
  v = 2.0 * Math.atan(Math.sqrt((1.0 + e) / (1.0 - e)) * Math.tan(E / 2.0));
  v = (PI2 + v + ap) % PI2; /// True Anomaly
  o = id.o;
  e = Math.sin(v);
  v = Math.cos(v);
  n = Math.sin(o);
  o = Math.cos(o);
  M = Math.cos(id.i);
  x[0] = r * (o * v - (n * e * M));
  x[1] = -r * e * Math.sin(id.i);
  x[2] = r * (n * v + (o * e * M));
  x = popmatrix(orbitmatrix,x);
  return x;

Here's the rest and the proof in 3 method comparison- plus Walker's table and Azimuth. http://innerbeing.epizy.com/charts/orbital.html Look no further.

if you just want Gravity:

function gravitate(i) {
  var p = [,] /// old pnt
  var v = [,] /// old velocity;
  var cxy = [cx,cy]; /// Center
  var dv = subv2(cxy,p);
  var f = Math.pow(dv[0],2) + Math.pow(dv[1],2);
  f = Mass / (f * Math.sqrt(f));
  cxy = multv2(dv,f); v = addv2(v,cxy);
  newvel = v;
  newpnt = addv2(p,v);

One more thing, if your planet moved, orbiters need to be updated first. addv(orbiterpoint,subv(new,old)) - gravitate() assumes stationary center.


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.