Laws of Planetary Motion

First several packages are required:

--> load(draw)
--> load(newton1)

An ellipse is determined by the length 2a of its major axis and the eccentricity e (0 < e < 1). The trajectory $(x,y) = (a\cos E,a\sqrt{1-e^2}\sin E)$ is drawn by using the parameter E, called the eccentric anomaly.

--> a: 2; e: 0.6;
--> DRAW: [grid=true, xrange=[-3,4], yrange=[-3,3],
       color=blue,
       parametric(a*cos(E),a*sqrt(1-e^2)*sin(E),E,0,2*%pi)];
    apply(draw2d, DRAW);

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

--> ff: [a*e,0]
--> DRAW: append(DRAW,
      [color=red,
       point_size=2, point_type=filled_circle, points([ff])]);
    apply(draw2d, 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 obtained, in which the mean anomaly M moves equally in time. Here we need to calculate the eccentric anomaly E from the equally separated values of M by solving the inverse Kepler's equation via Newton's method.

--> 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]),a*sqrt(1-e^2)*sin(E[i])]-ff,i,1,12)
--> DRAW: append(DRAW,
             append([color=red, head_length=0.05],
                    makelist(vector(ff,RR[i]),i,1,12)));
    apply(draw2d, DRAW);

Determining the velocity of the planet. Let $\mathbf{r}(t)$ be the trajectory of the planet viewed from the sun as the origin of coordinate. It satisfies the law of gravitation $\mathbf{r}''(t) = -GM(\mathbf{r}(t)/\vert\mathbf{r}(t\vert^3)$ , and the velocity $\mathbf{r}'(t)$ can be determined on the particular elliptical orbit. For the sake of simplicity we assume that GM = 1. From the conservation of angular momentum $\mathbf{r}(t)\times\mathbf{r}'(t) = \mathbf{h}$ with a constant vector $\mathbf{h}$ perpendicular to the orbit, we obtain the equation of velocity

$\displaystyle \mathbf{r}'(t)\times\mathbf{h} = \mathbf{u}(t) + \mathbf{c}$

where $\mathbf{u}(t) = \mathbf{r}(t)/\vert\mathbf{r}(t)\vert$ and $\mathbf{c}$ is another constant vector. Here we draw the direction $\mathbf{u}(t) + \mathbf{c}$ .

--> UU: makelist(rr/sqrt(rr.rr),rr,RR)
--> cc: [e,0]
--> VVH: makelist(uu+cc,uu,UU)
--> DRAW: append(DRAW,
              append([color=yellow, head_length=0.1],
                     makelist(vector(RR[i]+ff,VVH[i]),i,1,12)));
    apply(draw2d, DRAW)

It is deduced that $\mathbf{r}(t)$ satisfies the equation of trajectory

$\displaystyle a(1-e^2) = \vert\mathbf{h}\vert^2 = \mathbf{r}(t)\cdot(\mathbf{r}'(t)\times\mathbf{h}) = \mathbf{r}(t)\cdot(\mathbf{u}(t)+\mathbf{c})$

From the equation of velocity and of trajectory we can determine the direction of $\mathbf{r}'(t)$ and the magnitude $\vert\mathbf{r}'(t)\vert$ .

--> 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]+ff,VV[i]),i,1,12)));
    apply(draw2d, DRAW)

Here you can download the entire demonstration.

kepler.wxm