# 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? – Lars Viklund Dec 3 '15 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 – planetmaker Dec 3 '15 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. – imnegan Dec 3 '15 at 11:02
• Planetmaker- thanks for the recommendation, I'll check it out. – imnegan Dec 3 '15 at 11:03

## 1 Answer

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);


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? – imnegan Apr 17 at 12:17