# Orbital mechanics: orbit as a function of time. Universal variable formulation?

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.

Cheers!

• 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? Dec 3, 2015 at 8:27
• 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 Dec 3, 2015 at 10:02
• 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. Dec 3, 2015 at 11:02
• Planetmaker- thanks for the recommendation, I'll check it out. Dec 3, 2015 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);
else
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)));

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)
break;
}
return ret;
}


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 :

//Position
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;

//Velocity
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

Orbital elements

Kepler's equation

True anomaly

Eccentricity vector

Argument of periapsis

Newton's method for numerical approximation

• 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? Apr 17, 2019 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:

//2D
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;