From Dancing Planet to Kepler’s Laws

This most beautiful system of the sun, planets, and comets, could only proceed from the counsel and dominion of an intelligent powerful Being” Sir. Issac Newton

When I was seven years old, I had the notion that all planets dance around the sun along a wavy orbit (see Fig. 1).

Fig. 1

Many years later, I took on a challenge to show mathematically the orbit of my ‘dancing planet’ . This post is a long overdue report of my journey.

Shown in Fig. 2 is the sun and a planet in a x-y-z coordinate system. The sun is at the origin. The moving planet’s position is being described by x=x(t), y=y(t), z=z(t).

Fig. 2 r=\sqrt{x^2+y^2+z^2}, F=G\frac{Mm}{r^2}, F_z=-F\cos(c)=-F\cdot\frac{z}{r}

According to Newton’s theory, the gravitational force sun exerts on the planet is

F=-G\cdot M \cdot m \cdot \frac{1}{r^2}(\frac{x}{r},\frac{y}{r}, \frac{z}{r})=-\mu\cdot m \cdot \frac{1}{r^3}\cdot(x, y, z)

where G is the gravitational constant, M, m the mass of the sun and planet respectively. \mu = G\cdot M.

By Newton’s second law of motion,

\frac{d^2x}{dt^2} = -\mu\frac{x}{r^3},\quad\quad\quad(0-1)

\frac{d^2y}{dt^2} = -\mu\frac{y}{r^3},\quad\quad\quad(0-2)

\frac{d^2z}{dt^2} = -\mu\frac{z}{r^3}.\quad\quad\quad(0-3)

y \cdot (0-3) - z \cdot (0-2) yields

y\frac{d^2z}{dt^2}-z\frac{d^2y}{dt^2} = -\mu\frac{yz}{r^3}+ \mu\frac{yz}{r^3}=0.

Since

y\frac{d^2z}{dt^2}-z\frac{d^2y}{dt^2} = \frac{dy}{dt}\frac{dz}{dt}+y\frac{d^2z}{dt^2}-\frac{dz}{dt}\frac{dy}{dt}-z\frac{d^2y}{dt^2}=\frac{d}{dt}(y\frac{dz}{dt}-z\frac{dy}{dt}),

it must be true that

\frac{d}{dt}(y\frac{dz}{dt}-z\frac{dy}{dt}) = 0.

i.e.

y\frac{dz}{dt}-z\frac{dy}{dt}=A\quad\quad\quad(0-4)

where A is a constant.

Similarly,

z\frac{dx}{dt}-x\frac{dz}{dt}= B,\quad\quad\quad(0-5)

x\frac{dy}{dt}-y\frac{dx}{dt}= C\quad\quad\quad(0-6)

where B,C are constants.

Consequently,

Ax=xy\frac{dz}{dt} - xz\frac{dy}{dt},

By=yz\frac{dx}{dt} - xy\frac{dz}{dt},

Cz=xz\frac{dy}{dt}-yz\frac{dx}{dt}.

Hence

Ax + By +Cz=0.\quad\quad\quad(0-7)

If C \ne 0 then by the following well known theorem in Analytic Geometry:

If A, B, C and D are constants and A, B, and C are not all zero, then the graph of the equation Ax+By+Cz+D=0 is a plane“,

(0-7) represents a plane in the x-y-z coordinate system.

For C=0, we have

\frac{d}{dt}(\frac{y}{x})=\frac{x\frac{dy}{dt}-y\frac{dx}{dt}}{x^2}\overset{(0-6)}{=}\frac{C}{x^2}=\frac{0}{x^2}=0.

It means

\frac{y}{x}=k

where k is a constant. Simply put,

y=k x.

Hence, (0-7) still represents a plane in the x-y-z coordinate system (see Fig. 3(a)).

Fig. 3

The implication is that the planet moves around the sun on a plane (see Fig. 4).

Fig. 4

By rotating the axes so that the orbit of the planet is on the x-y plane where z \equiv 0 (see Fig. 3), we simplify the equations (0-1)-(0-3) to

\begin{cases} \frac{d^2x}{dt^2}=-\mu\frac{x}{r^3} \\ \frac{d^2y}{dt^2}=-\mu\frac{y}{r^3}\end{cases}.\quad\quad\quad(1-1)

It follows that

\frac{d}{dt}((\frac{dx}{dt})^2 + (\frac{dy}{dt})^2)

= 2\frac{dx}{dt}\cdot\frac{d^2x}{dt^2} + 2 \frac{dy}{dt}\cdot\frac{d^2y}{dt^2}

\overset{(1-1)}{=}2\frac{dx}{dt}\cdot(-\mu\frac{x}{r^3})+ 2\frac{dy}{dt}\cdot(-\mu\frac{y}{r^3})

= -\frac{\mu}{r^3}\cdot(2x\frac{dx}{dt}+2y\frac{dy}{dt})

=  -\frac{\mu}{r^3}\cdot\frac{d(x^2+y^2)}{dt}

= -\frac{\mu}{r^3}\cdot\frac{dr^2}{dt}

= -\frac{\mu}{r^3} \cdot 2r \cdot \frac{dr}{dt}

= -\frac{2\mu}{r^2} \cdot \frac{dr}{dt}.

i.e.,

\frac{d}{dt}((\frac{dx}{dt})^2 + (\frac{dy}{dt})^2) = -\frac{2\mu}{r^2} \cdot \frac{dr}{dt}.

Integrate with respect to t,

(\frac{dx}{dt})^2+(\frac{dy}{dt})^2 = \frac{2\mu}{r} + c_1\quad\quad\quad(1-2)

where c_1 is a constant.

We can also re-write (0-6) as

x\frac{dy}{dt}-y\frac{dx}{dt}=c_2\quad\quad\quad(1-3)

where c_2 is a constant.

Using polar coordinates

\begin{cases} x= r\cos(\theta) \\ y=r\sin(\theta) \end{cases},

Fig. 5

we obtain from (1-2) and (1-3) (see Fig. 5):

(\frac{dr}{dt})^2 + (r\frac{d\theta}{dt})^2-\frac{2\mu}{r} = c_1,\quad\quad\quad(1-4)

r^2\frac{d\theta}{dt} = c_2.\quad\quad\quad(1-5)

If the speed of planet at time t is v then from Fig. 6,

Fig. 6

v = \lim\limits_{\Delta t \rightarrow 0}\frac{\Delta l}{\Delta t} = \lim\limits_{\Delta t\rightarrow 0}\frac{l_r}{\Delta t}\overset{l_r=r\Delta \theta}{=}\lim\limits_{\Delta t \rightarrow 0}\frac{r\cdot \Delta \theta}{\Delta t}=r\cdot\lim\limits_{\Delta t\rightarrow 0}\frac{\Delta \theta}{\Delta t}=r\cdot\frac{d\theta}{dt}

gives

v = r\frac{d\theta}{dt}.\quad\quad\quad(1-6)

Suppose at t=0, the planet is at the greatest distance from the sun with r=r_0, \theta=0 and speed v_0. Then the fact that r attains maximum at t=0 implies (\frac{dr}{dt})_{t=0}=0. Therefore, by (1-4) and (1-5),

(\frac{dr}{dt})^2_{t=0} + (r\frac{d\theta}{dt})^2_{t=0}-\frac{s\mu}{r} = 0+ v_0^2-\frac{2\mu}{r}=c_1,

r (r\frac{d\theta}{dt})_{t=0}=r_0v_0=c_2.

i.e.,

c_1=v_0^2-\frac{2\mu}{r_0},\quad\quad\quad(1-7)

c_2=v_0 r_0.\quad\quad\quad(1-8)

We can now express (1-4) and (1-5) as:

\frac{dr}{dt} = \pm \sqrt{c_1+\frac{2\mu}{r}-\frac{c_2^2}{r^2}},\quad\quad\quad(1-9)

\frac{d\theta}{dt} = \frac{c_2}{r^2}.\quad\quad\quad(1-10)

Let

\rho = \frac{c_2}{r}\quad\quad\quad(1-11)

then

\frac{d\rho}{dr} = -\frac{c_2}{r^2},\quad\quad\quad(1-12)

r=\frac{c_2}{\rho}.\quad\quad\quad(1-13)

By chain rule,

\frac{d\theta}{dt} = \frac{d\theta}{d\rho}\cdot\frac{d\rho}{dr}\cdot\frac{dr}{dt}.

Thus,

\frac{d\theta}{d\rho} = \frac{\frac{d\theta}{dt}}{ \frac{d\rho}{dr} \cdot \frac{dr}{dt}}

\overset{(1-10), (1-12), (1-9)}{=} \frac{\frac{c_2}{r^2}}{ (-\frac{c_2}{r^2})\cdot(\pm\sqrt{c_1+\frac{2\mu}{r}-\frac{c_2^2}{r^2}}) }

\overset{(1-11)}{=} \mp\frac{1}{\sqrt{c_1-\rho^2+2\mu(\frac{\rho}{c_2})}}

= \mp\frac{1}{\sqrt{c_1+(\frac{\mu}{c_1})^2-\rho^2+2\mu(\frac{\rho}{c_2}) -(\frac{\mu}{c_2})^2}}

= \mp\frac{1}{\sqrt{c_1+(\frac{\mu}{c_1})^2-(\rho^2-2\mu(\frac{\rho}{c_2}) +(\frac{\mu}{c_2})^2)}}

= \mp\frac{1}{\sqrt{c_1+(\frac{\mu}{c_1})^2-(\rho-\frac{\mu}{c_2})^2}}.

That is,

\frac{d\theta}{d\rho} = \mp\frac{1}{\sqrt{c_1+(\frac{\mu}{c_1})^2-(\rho-\frac{\mu}{c_2})^2}}.\quad\quad\quad(1-14)

Since

c_1+(\frac{\mu}{c_2})^2\overset{(1-7)}{=}v_0^2-\frac{2\mu}{r_0}+(\frac{\mu}{v)0r_0})^2=(v_0-\frac{\mu}{v_0r_0})^2,

we let

\lambda = \sqrt{c_1 + (\frac{\mu}{c_2})^2}=\sqrt{(v_0-\frac{\mu}{v_0r_0})^2}=|v_0-\frac{\mu}{v_0r_0}|.

Notice that \lambda \ge 0.

By doing so, (1-14) can be expressed as

\frac{d\theta}{d\rho} =\mp \frac{1}{\sqrt{\lambda^2-(\rho-\frac{\mu}{c_2})^2}} .

Take the first case,

\frac{d\theta}{d\rho} = -\frac{1}{\sqrt{\lambda^2-(\rho-\frac{\mu}{c_2})^2}} .

Integrate it with respect to \rho gives

\theta + c = \arccos(\frac{\rho-\frac{\mu}{c_2}}{\lambda})

where c is a constant.

When r=r_0, \theta=0,

c = \arccos(1)=0 or c = \arccos(-1) = \pi.

And so,

\theta = \arccos(\frac{\rho-\frac{\mu}{c_2}}{\lambda}) or \theta+\pi = \arccos(\frac{\rho-\frac{\mu}{c_2}}{\lambda}).

For c = 0,

\lambda\cos(\theta) = \rho-\frac{\mu}{c_2}.

By (1-11), it is

\frac{c_2}{r}-\frac{\mu}{c_2} = \lambda \cos(\theta).\quad\quad\quad(1-15)

Fig. 7

Solving (1-15) for r yields

r=\frac{c_2^2}{c_2 \lambda \cos(\theta)+\mu}=\frac{\frac{c_2^2}{\mu}}{\frac{c_2}{\mu}\lambda \cos(\theta)+1}\overset{p=\frac{c_2^2}{\mu}, e=\frac{c_2 \lambda}{\mu}}{=}\frac{p}{e \cos(\theta) + 1}.

i.e.,

r = \frac{p}{e \cos(\theta) + 1}.\quad\quad\quad(1-16)

Studies in Analytic Geometry show that for an orbit expressed by (1-16), there are four cases to consider depend on the value of e:

We can rule out parabolic and hyperbolic orbit immediately for they are not periodic. Given the fact that a circle is a special case of an ellipse, it is fair to say:

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

In fact, this is what Kepler stated as his first law of planetary motion.

Fig. 8

For c=\pi,

\theta + \pi = \arccos(\frac{\rho-\frac{\mu}{c_2}}{\lambda})

from which we obtain

r=\frac{c_2^2}{c_2 \lambda \cos(\theta+\pi)+\mu}=\frac{\frac{c_2^2}{\mu}}{\frac{c_2}{\mu}\lambda\cos(\theta+\pi)+1}\overset{p=\frac{c_2^2}{\mu}, e=\frac{c_2 \lambda}{\mu}}{=}\frac{p}{e \cos(\theta+\pi) + 1}.\quad\quad(1-17)

This is an ellipse. Namely, the result of rotating (1-16) by hundred eighty degrees or assuming r attains its minimum at t=0.

The second case

\frac{d\theta}{d\rho} = +\frac{1}{\sqrt{\lambda^2-(\rho-\frac{\mu}{c_2})^2}}

can be written as

-\frac{d\theta}{d\rho} = -\frac{1}{\sqrt{\lambda^2-(\rho-\frac{\mu}{c_2})^2}} .

Integrate it with respect to \rho yields

-\theta + c = \arccos(\frac{\rho-\frac{\mu}{c_2}}{\lambda})

from which we can obtain (1-16) and (1-17) again.

Fig. 9

Over the time duration \Delta t, the area a line joining the sun and a planet sweeps an area A (see Fig. 9).

A = \int\limits_{t}^{t+\Delta t}\frac{1}{2}r\cdot v\;dt \overset{(1-6)}{=} \int\limits_{t}^{t+\Delta t}\frac{1}{2}r\cdot r\frac{d\theta}{dt}\;dt=\int\limits_{t}^{t+\Delta t}\frac{1}{2}r^2\frac{d\theta}{dt}\;dt\overset{(1-5)}{=}\int\limits_{t}^{t+\Delta t}\frac{1}{2}c_2\;dt = \frac{1}{2}c_2\Delta t.

It means

A = \frac{1}{2}c_2\Delta t\quad\quad\quad(2-1)

or that

\frac{A}{\Delta t} = \frac{1}{2}c_2

is a constant. Therefore,

A line joining the Sun and a planet sweeps out equal areas during equal intervals of time.

This is Kepler’s second law. It suggests that the speed of the planet increases as it nears the sun and decreases as it recedes from the sun (see Fig. 10).

Fig. 10

Furthermore, over the interval T, the period of the planet’s revolution around the sun, the line joining the sun and the planet sweeps the enire interior of the planet’s elliptical orbit with semi-major axis a and semi-minor axis b. Since the area enlosed by such orbit is \pi ab (see “Evaluate a Definite Integral without FTC“), setting \Delta t in (2-1) to T gives

{\frac{1}{2}c_2 T} = {\pi a b} \implies  {\frac{1}{4}c_2^2 T^2}={\pi^2 a^2 b^2} \implies T^2 = \frac{4\pi^2 a^2 b^2}{c_2^2} \implies \frac{T^2}{a^3} = \frac{4\pi^2b^2}{c_2^2a}. (3-1)

While we have p = \frac{c_2^2}{\mu} in (1-16), it is also true that for such ellipse, p=\frac{b^2}{a} (see “An Ellipse in Its Polar Form“). Hence,

\frac{b^2}{a}=\frac{c_2^2}{\mu}\implies c_2^2=\frac{\mu b^2}{a}.\quad\quad\quad(3-2)

Substituting (3-2) for c_2^2 in (3-1),

\frac{T^2}{a^3} = \frac{4\pi^2 b^2}{(\frac{\mu b^2}{a})a}=\frac{4\pi^2}{\mu} \overset{\mu=GM}{=}\frac{4\pi^2}{GM}.\quad\quad\quad(3-3)

Thus emerges Kepler’s third law of planetary motion:

The square of the orbital period of a planet is directly proportional to the cube of the semi-major axis of its orbit.

Established by (3-3) is the fact that the proportionality constant is the same for all planets orbiting around the sun.

This image has an empty alt attribute; its file name is main-qimg-e7446d7decf4bdc83f8b712109470213.jpg

1 thought on “From Dancing Planet to Kepler’s Laws

  1. Pingback: An Ellipse in Its Polar Form | Vroom

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s