What are the equations for an orbit with time as parameter?

Question

I am preparing a talk on the Equation of Time, and I want to produce an animation of the movement of the earth around the sun, in an exaggerated ellipse, to illustrate the speeding up and slowing down in different parts of the orbit. Ideally I would like to have x and y as functions of time, e.g x = f(t) and y = g(t). Accuracy is not the prime consideration, a reasonable approximation will do.

I will use a Processing sketch to produce the animation.

Can anyone please tell me what these equations are?

Answer 1

I have done a lot more research into this and found an answer. It comes down to solving Kepler's Equation, and my main references are: Kepler's equation and the Equation of Centre which gives an iterative method for solving the equation, and Determination of Position in Orbits which has the formulae for calculating orbital coordinates.

Here is the Processing sketch. If you are not familiar with Processing, it is a language developed for graphic designers, but much more versatile than just that, based on Java and hence C++. It expects two functions, one is setup(), called once at the start up of the sketch, and a second draw() that is repeated over and over until the sketch is stopped. In case you wonder: the whole is encased in a main() before it is compiled. For more details see the Processing home page.

int wsize = 400;
float rp = 9.0; // radius of planet symbol
float rs = 15.0; // --"-- sun
// Constants of orbital ellipse 
float a = 150, b; // major and minor axes
float e = 0.5; // eccentricity
float halfa = a/2.0;
float c; // focal distance from centre
float E, E0; // eccentric anomaly, temporary values in iteration
float v; // true anomaly = polar coord angle
float M; /* mean anomaly which is effectively the ..
            .. time variable round the orbit*/
float r, x, y; // polar coord radius, rect coords
float k, cosE, cosv; // useful intermediate variables
float twopi = 2*PI;
float PIby2 = PI/2.0;
float Minc = twopi/250; // time proxy increment

void setup()
{ size(wsize,wsize);
  //e = sqrt(((a*a)-(b*b))/(a*a)); kept for reference
  b = a*sqrt(1-e*e);
  c = e*a;
  k = 1-e*e;
  M = 0; // Mean anomaly increases linearly with time
}

void draw()
{ background(255);
  translate(wsize/2,wsize/2); // move origin to centre of display window
  ellipse(0,0,2*a,2*b); // show orbit
  spot(c,0,rs,#FFF303); // draw sun
  translate(c,0); // move origin to focus
  // Calculate Eccentic Anomaly by iteration
  E = M; E0 = 1000; // make sure while starts
  while(abs(E-E0) > 0.000001)
  { E0 = M + e*sin(E);
    E = E0;
  } 
  // calculate polar coordinates of planet
  cosE= cos(E);
  cosv = (cosE-e)/(1-e*cosE);
  v= acos(cosv);
  r = a*(1-e*e)/(1+e*cosv);
  // calculate x,y coordinates of planet
  x = r * cosv;
  y = -r*sin(v);// negative makes it go anti-clockwise
  if (M > PI){y = -y;} // show second half of orbit
  // println(v+" : "+x+" : "+y);
  // line(0,0,x,y);
  spot(x,y, rp, #03ADFF); // draw planet
  M += Minc;
  if(M > twopi){M = 0.0;}
}

void spot(float x, float y, float w, color tint)
{ fill(tint);
  ellipse(x,y,w,w);
  fill(255); // restore white background
}