Kepler's laws of planetary motion

load(draw);

load(newton1);

An ellipse is determined by the length 2a of its major axis and
the eccentricity e (0 < e < 1).

a: 2; e: 0.6;

In polar coordinates the ellipse is determined by (r, theta)
satisfying the polar equation:

r: a*(1-e^2)/(1+e*cos(theta));

DRAW: [grid=true, xrange=[-5,2], yrange=[-3,3],
       color=blue,
       parametric(r*cos(theta),r*sin(theta),theta,0,2*%pi)];
apply(wxdraw2d, DRAW);

Kepler's first law: A planet revolves around the sun in an elliptical orbit
with the sun at the focus (0,0).

ff: [0,0];

DRAW: append(DRAW,
      [color=orange,
       point_size=2, point_type=filled_circle, points([ff])]);
apply(wxdraw2d, DRAW);

Kepler's second law: The line joining the sun to a planet sweeps out equal areas
in equal times. Furthermore, the Kepler's equation M = E - e sin(E) is observed,
in which the mean anomaly M moves equally in time.

for i from 1 thru 12 do (E[i]: newton(x-e*sin(x)-3.14*(i/6), x, 3.14*(i/6), 1/50));

RR: makelist([a*(cos(E[i])-e),a*sqrt(1-e^2)*sin(E[i])],i,1,12);

DRAW: append(DRAW,
         append([color=orange, head_length=0.05],
                makelist(vector(ff,RR[i]),i,1,12)));
apply(wxdraw2d, DRAW);

Assuming GM = 1, we obtain cc = [e,0]. Then obtain the direction orthogonal
to the velocity of the planet.

UU: makelist(rr/sqrt(rr.rr),rr,RR);

cc: [e,0];

VVH: makelist(uu+cc,uu,UU);

DRAW: append(DRAW,
         append([color=red, head_length=0.1],
                makelist(vector(RR[i],VVH[i]),i,1,12)));
apply(wxdraw2d, DRAW);

Calculate the velocity of the planet. The magnitude |r'| of velocity r' is
determined by |r'|*|h| = |u + c| with |h| = sqrt(a*(1-e^2)).

VV: makelist([-dd[2],dd[1]]/sqrt(a*(1-e^2)), dd, VVH);

DRAW: append(DRAW,
         append([color=green, head_length=0.1],
                makelist(vector(RR[i],VV[i]),i,1,12)));
apply(wxdraw2d, DRAW);