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
}
No comments:
Post a Comment