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.


You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .