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

Question

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!

Answer 1

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

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

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

Answer 2

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