Kepler's Laws

AP1603
M. E. Mauel

Kepler

The German astronomer Johannes Kepler (b.  Dec.  27, 1571, d.  Nov.  15, 1630) was the first strong supporter of Copernicus's view that the planets revolved the Sun and the discoverer of the famous three laws of planetary motion.  Kepler had planned to enter religious life, attended seminaries at Adelberg and Maulbronn, and studied theology, philosophy, and mathematics at the University of Tubingen.  At Tubingen, Kepler's scientific ability attracted notice. In his early 20's, Kepler became a professor of mathematics and astronomy at Graz.
At the age of 24, Kepler published Mysterium cosmographicum (Cosmographic Mystery, 1596), in which he defended the Copernican theory and described his ideas on the structure of the planetary system.  Influenced by the Pythagoreans, Kepler viewed the universe as being governed by geometric relationships that conform to the inscribed and circumscribed circles of the five regular polygons.
Although he was not a Copernican himself, Tycho Brahe, the mathematician at the court of Emperor Rudolph II at Prague, was so impressed with Kepler's work that in 1600 he invited Kepler to come to Prague as his assistant.  Confronted with the Catholic persecution of the Protestant minority in Graz, Kepler gladly accepted.  When Brahe died the following year, Kepler was appointed his successor and thus inherited Brahe's scientific legacy.
This legacy included many accurate positional determinations of the planets, especially those of Mars.  Kepler now embarked on an intensive study of the true orbits of the planets.  Abandoning the ancient belief that the planets must move in perfect circles, Kepler concentrated on Mars.  He proved that the orbit of Mars is an ellipse, with the Sun occupying one of its two foci.  This, the first of Kepler's laws of planetary motion, appeared in Astronomia nova (New Astronomy) in 1609, with the second "law of areas" governing planetary velocity.
Always guided by the concept of beauty in the structure of the universe, and specifically by a theory of harmony in geometric figures, numbers, and music, Kepler, in his Harmonices mundi (Harmonies of the World, 1619), announced his third law--a relationship between the orbital periods and the distances of the planets from the Sun.  His belief that the Sun regulates the velocity of the planets was a milestone in scientific thought, laying the foundation for Newton's theory of universal gravitation.
Among Kepler's numerous scientific contributions are an influential treatise on the theory of optics (1604), a treatise on optics as applied to telescope lenses (1611), a work offering physical explanations of the appearance of a nova in 1604, and an enthusiastic acceptance of and elaboration on Galileo's observations with a telescope (1610).  His Epitome astronomiae Copernicanae (Introduction to Copernican Astronomy, 1618-21) became one of the most widely read treatises on astronomy in Europe.  Kepler's last great work, known as the Rudolphine Tables (1627), was a widely used compilation of accurate tables of planetary motion.
The posthumous Somnium (Dream, 1634), on which Kepler labored until shortly before his death, is indicative of his fertile mind.  In this work, Kepler describes a journey to the Moon and discusses the existence of lunar inhabitants.  
(From Steven J. Dick, Groiler's Encyclopedia.)

Kepler's Three Laws

Kepler's three laws of planetary motion describe the shape of planetary motion on the plane, the velocities of the planets, and the distances of the planets from the Sun. The laws are:
1. All planets move in elliptic orbits with the Sun at one focus.
2. The line joining a planet to the Sun sweeps out equal areas in equal times.
3. If T is equal to the orbital period of a planet (i.e. the time required to complete one orbit), and if a is the semimajor axis of the orbit (i.e. one half the long distance across the ellispe), then the ratio (T^2)/(a^3) is a constant.

Planetary Motion about a Massive Star

According to Sir Issac Newton, the force between a planet and the Sun is
	(Subscript[F, G]) = G (Subscript[M, s])m / (r^2)
where m is the mass of the planet and (r^2)= (x^2)+ (y^2) is the square of the distance between the planet and the Sun. 
The force of gravity always attracts one object towards the other. Force is a vector, and as the planet orbits the Sun, the direction of the force also revolves. From the diagram above, the x-directed force is equal to the total force, (Subscript[F, G]), times the ratio, x/r. The y-component of the force is y/r times (Subscript[F, G]).
Newton also explained the force on the Sun due to the planet is equal and opposite to the force on the planet from the Sun. However, if the Sun is much, much more massive than the planet, the Sun does not significantly move as compared to the planet. In this approximation, we need only follow the trajectory of the planet. (Physics student's should note: Newton also showed how to describe the motion of any two objects interacting through gravity, like binary stars, and this description is often called Kepler's ammended laws.)
In the x-y plane, the equations of motion of a planet about a massive star are:
	x'[t] = vx[t]
	y'[t] = vy[t]
	vx'[t] = - G(Subscript[M, S]) x / (r^3)
	vy'[t] = - G(Subscript[M, S]) y / (r^3)
where the prime denotes teh first derivative with respect to time, df /dt = f'.
The product of Newton's gravitational constant and the mass of the sun, G(Subscript[M, S]), takes on a simple form if we measure distance in astronomical units, or AU, and if we measure time in years. In these units, G(Subscript[M, S]) = 4(π^2) (AU^3)/(yr^2). One AU is equal to the average distance between the Sun and the Earth, about 1.5 × (10^11)m.

Equations of Motion

The equations of motion for a planet describe the time evolution of the planet's phase-space location and its initial position.

p = {x[t], y[t], vx[t], vy[t]};

dpdt = {D[x[t],t] == vx[t],
D[y[t],t] == vy[t],
D[vx[t],t] == -4 (π^2) x[t]/r^3,
D[vy[t],t] == -4 (π^2) y[t]/r^3} /. r -> Sqrt[x[t]^2 + y[t]^2]

p0 = {x[0] == x0, y[0] == y0,vx[0] == vx0, vy[0] == vy0}

eqs = dpdt ~ Join ~p0

Earth's Orbit

In astronomical units (AUand year), the trajectory of the Earth is relatively easy to calculate. The key is finding the appropriate initial condition. 
Let's set the position of Earth to be {x, y} = {1, 0}. The initial velocity is approximately the circumference, 2πr = 2π, dividedby 1 year, or {vx, vy} = {0, 2π}.

In Cartesian coordinates, the equations of motion are two difficult for Mathematica to solve using DSolve. (If we do some algebra and transformation of variables, we can use DSolve for the two-body central-force problem, but we'll discuss this later.)

earth = DSolve[eqs /.x0 -> 1 /. y0 -> 0 /. vy0->2 π/. vx0 -> 0, p, t]

Instead, we must use NDSolve.

earth = NDSolve[eqs /.x0 -> 1 /. y0 -> 0 /. vy0->2 π/. vx0 -> 0, p, {t, 0, 1.2}]

ParametricPlot[Evaluate[{x[t],y[t]}/.earth],{t,0,1},AspectRatio->Automatic]

The orbit is a circle because of our choice of the initial condition.

Off[General::spell, General::spell1];

earthA = NDSolve[eqs /.x0 -> 1 /. y0 -> 0 /. vy0->2 π/1.1/. vx0 -> 0, p, {t, 0, 1.2}];
earthB = NDSolve[eqs /.x0 -> 1 /. y0 -> 0 /. vy0->2 π/1.2/. vx0 -> 0, p, {t, 0, 1.2}];
earthC = NDSolve[eqs /.x0 -> 1 /. y0 -> 0 /. vy0->2 π 1.1/. vx0 -> 0, p, {t, 0, 1.2}];
earthD = NDSolve[eqs /.x0 -> 1 /. y0 -> 0 /. vy0->2 π 1.2/. vx0 -> 0, p, {t, 0, 1.2}];

ParametricPlot[Evaluate[{{x[t],y[t]}/.earth,{x[t],y[t]}/.earthA,{x[t],y[t]}/.earthB,{x[t],y[t]}/.earthC,{x[t],y[t]}/.earthD}],{t,0,1},AspectRatio->Automatic,PlotStyle->{{},{Hue[(0/4)]},{Hue[(1/4)]},{Hue[(2/4)]},{Hue[(3/4)]}},Compiled->False]

Halley's Comet

Halley's coment has anorbital period of 76 yr and a closest distance tothe Sun of 0.59 AU. 
The challenge is finding the initial condition which determines the orbit.
Below, we try the "brute force" or "rial-and-error" approach--perfect for computers...

halley = NDSolve[eqs /.x0 -> 0.59 /. y0 -> 0 /. vy0->11.46/. vx0 -> 0, p, {t, 0, 76}] // Flatten

ParametricPlot[Evaluate[{x[t],y[t]}/.halley],{t,0,76}]

Show[%,AspectRatio->Automatic]

Plot[Evaluate[{vx[t],vy[t]}/.halley],{t,0,76},PlotStyle->{Red,Green}]

Notice that the orbital period is not quite 76 years. We'll try forming a table of possible orbits.

halleyTable = Table[NDSolve[eqs /.x0 -> 0.59 /. y0 -> 0 /. vy0->v/. vx0 -> 0, p, {t, 0, 76}] // Flatten, {v, 11.46, 11.48, 0.004}];

ParametricPlot[Evaluate[{x[t],y[t]}/.halleyTable],{t,36,40},PlotStyle->Table[{Hue[(i/5)]},{i,0,4}]]

So the  velocity of Halley's comet at its closest distance to the Sun is about v ≈ 11.472 AU/yr.

The Ellipse (optional)

(This section is optional. Students interested in techniques to use Mathematica to manipulate algebraic equations can look at the expressions below.)

An ellipse is is a curve such that the sum of the distances between each point on the ellipse to two fixed points (i.e. the foci) are constant.

A Centered, Symmetric Ellipse

If the two foci are located at {x, y} = {±c, 0}, and if the distance between the ellipse when it crosses the x axes is a, then the equation for the ellipse is

ellipse1 = Sqrt[(((x-c))^2)+ (y^2)] + Sqrt[(((x+c))^2)+ (y^2)] == 2 a

This can be turned into a more familiar form of the ellipse by replacing (c^2) with the difference between the major and minor axes.

eSol1 = Solve[ellipse1,y]//. (c^2)->(a^2) - b^2 // Simplify

ellipse2 = (x^2)/(a^2) + (y^2)/(b^2) == 1

eSol2 = Solve[ellipse2, y] // Simplify

Expand[y^2 /. eSol1[[1]]] == Expand[y^2/.eSol2[[1]]] (* QED *)

Plot[Evaluate[y/.eSol1/.b->(1/2)/.a->1],{x,-1,1},AspectRatio->Automatic]

The eccentricity, e = c / a, is defined to be the ratio of c to the semimajor axes, a. (This assumes the ordering, a > b and 0 ≥ c ≥ a.) From the replacements demonstrated above, the eccentricity is also given by the equation (b^2)= (a^2) (1 - (e^2)), or (a^2)(e^2)= (a^2) - (b^2).When e = 1 (and c = a), the ellipse becomes a line. When e= c = 0, the ellipse becomes a circle.

One Foci at the Origin

The orbit of a planet forms an ellipse with the Sun at one of the foci.

ellipse2 = Sqrt[(x^2)+ (y^2)] + Sqrt[(((x+2c))^2)+ (y^2)] == 2 a

eSol3 = Solve[ellipse2,y]//. c->Sqrt[(a^2) - b^2] // Simplify

c1 = Sqrt[(1^2)-(((1/2))^2)]  (* An example for a = 1 and b = 1/2 *)

N[c1]      (* c = e when a = 1 *)

{c1, 2c1, -1 -c1} // N

Plot[Evaluate[y/.eSol3/.b->(1/2)/.a->1],{x,-1-c1,1-c1},AspectRatio->Automatic]

For an ellipse, (y^2) varies as a quadratic in x

For an ellipse, the variation of (y^2)is given by a quadratic polynomial in x.

y^2 /. eSol3 // Collect[#,x]&

Kepler's First Law

Since we have the numerical solution to Newton's equations of motion, a relatively easy method to prove the solution is represented by an ellipse is to perform a least-squares fit between our (approximate numerical) solutions and a perfect ellipse.
For an ellipse, the variation of (y^2)is given by a quadratic polynomial in x. (See last section.)

Making a list of the coordinates of Halley once every year...

halleyEllipse = Table[{x[t],y[t]^2} /. halley, {t, 0, 76/2}];

ListPlot[halleyEllipse,PlotStyle->PointSize[0.02],AxesLabel->{"x","y^2"}]

The best fit ellipse to the numerical solution is

yfit[x_] = Sqrt[Fit[halleyEllipse, {1,x,x^2}, x]]

ParametricPlot[Evaluate[{{x[t],y[t]},{x[t],yfit[x[t]]}}/.halley],{t,0,(76/2)},PlotStyle->{{},{Red}}]

The error between the perfect elliptical fit and the numerical solution is pretty good:

Plot[Evaluate[y[t]-yfit[x[t]]/.halley],{t,0,(76/2)},PlotLabel->"Error (AU)"]

Kepler's Second Law

Kepler's second law states: The line joining a planet to the Sun sweeps out equal areas in equal times.

p1 = ParametricPlot[Evaluate[{x[t],y[t]} /. halley], {t, 0, 76}];

p2 = ParametricPlot[Evaluate[{x[t],y[t]} /. halley], {t,0, 2},PlotStyle->{Blue,Thickness[0.02]}];

p3 = ParametricPlot[Evaluate[{x[t],y[t]} /. halley], {t,35, 37},PlotStyle->{Green,Thickness[0.02]}];

Show[p1,p2,p3,AspectRatio->Automatic]

Kepler's second law is a statement of conservation of angular momentum.
The angular momentum about the Sun (located at the origin) is equal to 
	L = r × v
	L = x (Subscript[v, y]) - y (Subscript[v, x])
Expressed in this way (without the mass of the planet), L is equal to the incremental increase of the area spanned by the line connecting the planet to the star per incremental time.

l[t_] = x[t] vy[t] - y[t] vx[t] /. halley

Plot[l[t],{t,0,76}]

Plot[((100 (l[t]-l[0]))/(l[0])),{t,0,76},PlotLabel->"Percent Error of Angular Momentum"]

The total area of the ellipse is given by the integral of L around one complete orbit, A = πab.

Kepler's Third Law

Kepler's third law states:  If T is equal to the orbital period of a planet (i.e. the time required to complete one orbit), and if a is the semimajor axis of the orbit(i.e. one half the long distance across the ellispe), then the ratio (T^2)/(a^3) is a constant.

To prove Kepler's third law, we must develop a way to find the period of a planet's orbit and the length of its major axis.

Harmonic Constant for Halley

Plot[Evaluate[vx[t]/.halley],{t,0,76}]

Plot[Evaluate[vx[t]/.halley],{t,25,35}]

tHalley = 2t /. FindRoot[Evaluate[vx[t]/.halley], {t, 30}];

aHalley =(( x[t]/. halley /.t->0)- (x[t]/. halley/.t->tHalley/2))/2 ;

Harmonic Constant for Earth

Plot[Evaluate[vx[t]/.earth],{t,0,1}]

Plot[Evaluate[vx[t]/.earth],{t,0.8,1}]

tEarth = 1.0; aEarth = 1.0;

{tHalley^2/aHalley^3, tEarth^2/aEarth^3} // N

Summary

We've used NDSolve[…] to numerically integrate the equations of motion for an object orbiting a massive star. These solutions depend sensitively on the initial conditions of the object. We examined in detail the orbits of the planet Earth and the comet Halley. We showed that the solutions produced by Mathematica are in agreement with the three Laws of Kepler.