You played with an animation page of Satellite's orbit simulation (Kepler's laws). Now let's learn how these phenomena take place under the laws of physics.
Satellite's orbit is controlled by a gravity force. We treat the satellite as a particle without size and shape. We are showing an important fact (theorem):
If the gravity force comes from a mass with a perfect spherical shape, the orbit is mathematically solved and expressed by a quadratic curve (二次曲線) on a plane.In fact, all the orbit lines appearing on the animation was quadratic curves on a plane. It is the simplest model of a satellite's orbit. The goal of this page is to derive this from the law of physics.
Note that this theory is about a model (approximation) which is not the real orbit. The real orbit is not a true quadratic curve or on a true plane, because the satellite gets more complicated force than the gravity force coming from just a single source with a perfect sphere.We start with Newton's equation: \begin{eqnarray} {\bf F}=m {\bf a}\label{eq:Newton01} \end{eqnarray} First, we rewrite \eqref{eq:Newton01} in to a form using a plane polar coordinate \(r, \theta\). Namely, any vectors are represented by a linear combination of two unit vectors: \begin{eqnarray} {\bf e}_{r} = (\cos \theta, \sin \theta)\label{eq:er_cossin} \end{eqnarray} which directs from the origin to the point we are considering, and \begin{eqnarray} {\bf e}_{\theta} = (-\sin \theta, \cos \theta) \end{eqnarray} which directs perpendicular to \({\bf e}_{r}\) and anti-clockwise direction in a rotational movement whose center is at the origin. Remember that \({\bf e}_{r}\) and \({\bf e}_{\theta}\) change in response to the change of the position we consider, namely, the position of the satellite.
By definition, it is obvious that the position vector of the satellite \({\bf r}\) is described by \begin{eqnarray} {\bf r}=r{\bf e}_{r}\label{eq:def_r} \end{eqnarray} Then the velocity of the satellite \({\bf v}\) is: \begin{eqnarray} {\bf v}:=\dot{\bf r}\label{eq:def_v} \end{eqnarray} Note that the dot at the head of a letter means derivative with respect to time, nemaly, \(d/dt\).
By using \eqref{eq:def_r} and \eqref{eq:def_v}, we get more specific form of \({\bf v}\): \begin{eqnarray} {\bf v}:=\frac{d{\bf r}}{dt}=\frac{d(r{\bf e}_r)}{dt}=\dot{r}{\bf e}_r+r\dot{\bf e}_r\label{eq:calc_v01} \end{eqnarray} Because we have decided to describe all vectors with \({\bf e}_r\) and \({\bf e}_\theta\), we must describe \(\dot{\bf e}_r\) by them , too. Let's try using \eqref{eq:er_cossin}: \begin{eqnarray} \dot{\bf e}_r:&=&\frac{d{\bf e}_r}{dt}=\frac{d{\theta}}{dt}\frac{d{\bf e}_r}{d\theta}\nonumber\\ &=&\dot{\theta}\frac{d}{d\theta}(\cos \theta, \sin \theta)=\dot{\theta}(-\sin \theta, \cos\theta)\nonumber\\ &=&\dot{\theta}{\bf e}_\theta\label{eq:calc_dot_er04} \end{eqnarray} By substituting it to \eqref{eq:calc_v01}, we obtain: \begin{eqnarray} {\bf v}=\dot{r}{\bf e}_r+r\dot\theta{\bf e}_\theta\label{eq:calc_v03} \end{eqnarray} Then we proceed to satellite's acceleration \({\bf a}\): \begin{eqnarray} {\bf a}:&=&\dot{\bf v}=\frac{d}{dt}(\dot{r}{\bf e}_r+r\dot\theta{\bf e}_\theta)\nonumber\\ &=&(\ddot{r}{\bf e}_r+\dot{r}\dot{\bf e}_r+\dot{r}\dot\theta{\bf e}_\theta+r\ddot\theta{\bf e}_\theta+r\dot\theta\dot{\bf e}_\theta)\nonumber\\ &=&(\ddot{r}{\bf e}_r+\dot{r}\dot{\theta}{\bf e}_\theta+\dot{r}\dot\theta{\bf e}_\theta+r\ddot\theta{\bf e}_\theta+r\dot\theta\dot{\bf e}_\theta)\nonumber\\ &=&(\ddot{r}{\bf e}_r+2\dot{r}\dot{\theta}{\bf e}_\theta+r\ddot\theta{\bf e}_\theta+r\dot\theta\dot{\bf e}_\theta)\label{eq:calc_a01} \end{eqnarray} Here we can describe \(\dot{\bf e}_\theta\) by: \begin{eqnarray} \dot{\bf e}_\theta:&=&\frac{d{\bf e}_\theta}{dt}=\frac{d{\theta}}{dt}\frac{d{\bf e}_\theta}{d\theta}\nonumber\\ &=&\dot{\theta}\frac{d}{d\theta}(-\sin \theta, \cos \theta)=\dot{\theta}(-\cos \theta, -\sin\theta) \nonumber\\ &=&-\dot{\theta}{\bf e}_r\label{eq:calc_dot_etheta04} \end{eqnarray} By substituting it to \eqref{eq:calc_a01}, we obtain: \begin{eqnarray} {\bf a} &=&(\ddot{r}{\bf e}_r+2\dot{r}\dot{\theta}{\bf e}_\theta+r\ddot\theta{\bf e}_\theta-r\dot\theta^2{\bf e}_r)\nonumber\\ &=&(\ddot{r}-r\dot\theta^2){\bf e}_r+(2\dot{r}\dot{\theta}+r\ddot\theta){\bf e}_\theta\label{eq:calc_a03} \end{eqnarray} Meanwhile, by gravity's law, \begin{eqnarray} {\bf F}=-\frac{G M m}{r^{2}}{\bf e}_{r}\label{eq:gravity_law01} \end{eqnarray} By using \eqref{eq:Newton01}, \eqref{eq:calc_a03}, and \eqref{eq:gravity_law01}, we obtain: \begin{eqnarray} (\ddot{r}-r\dot\theta^2){\bf e}_r+(2\dot{r}\dot{\theta}+r\ddot\theta){\bf e}_\theta=-\frac{G M}{r^{2}}{\bf e}_{r} \end{eqnarray} Because \({\bf e}_{r}\) and \({\bf e}_\theta\) are linearly independent, \begin{eqnarray} &&\ddot{r}-r\dot\theta^2+\frac{G M}{r^{2}}=0\label{eq:dynamicequ04}\\ &&2\dot{r}\dot{\theta}+r\ddot\theta=0\label{eq:dynamicequ06} \end{eqnarray} By multiplying \(r\) to \eqref{eq:dynamicequ06}, \begin{eqnarray} 2r\dot{r}\dot{\theta}+r^2\ddot\theta=0\label{eq:dynamicequ08} \end{eqnarray} The left side of \eqref{eq:dynamicequ08} is equal to \(\frac{d}{dt}(r^2\dot\theta)\). Therefore, \begin{eqnarray} \frac{d}{dt}(r^2\dot\theta)=0\label{eq:dynamicequ10} \end{eqnarray} It means \(r^2\dot\theta\) is a constant (independent of time \(t\)). We represent it as a constant \(\alpha\). Namely, \begin{eqnarray} \alpha:=r^2\dot\theta\label{eq:def_alpha} \end{eqnarray} or \begin{eqnarray} \dot\theta = \frac{\alpha}{r^2}\label{eq:dynamicequ12} \end{eqnarray} By substituting \(\dot\theta\) in \eqref{eq:dynamicequ04} with \eqref{eq:dynamicequ12}, \begin{eqnarray} \ddot{r}-\frac{\alpha^{2}}{r^{3}}+\frac{G M}{r^{2}}\label{eq:dynamicequ14}=0 \end{eqnarray} Now we are going to rewrite this equation by changing the variable from \(t\) to \(\theta\). First, we rewrite \(\dot r\) and \(\ddot r\) as follows: \begin{eqnarray} \dot{r}&=&\frac{d r}{d t}=\frac{d r}{d \theta} \dot{\theta} = \frac{d r}{d \theta} \frac{\alpha}{r^{2}}\\ \ddot{r}&=&\frac{d \dot{r}}{d \theta} \dot{\theta}=\frac{d}{d \theta}\left(\frac{d r}{d \theta} \frac{\alpha}{r^{2}}\right) \dot{\theta}=\frac{\alpha}{r^{2}} \frac{d}{d \theta}\left(\frac{\alpha}{r^{2}} \frac{d r}{d \theta}\right)\label{\eq:ddotr_equation04} \end{eqnarray} (We used \eqref{eq:dynamicequ12} at the final step of both reforms.) Now we define \begin{eqnarray} u:=\frac{1}{r} \end{eqnarray} Then \begin{eqnarray} \frac{d u}{d \theta}=\frac{d}{d \theta} \frac{1}{r}=-\frac{1}{r^{2}} \frac{d r}{d \theta} \end{eqnarray} Therefore, \eqref{\eq:ddotr_equation04} becomes \begin{eqnarray} \ddot{r}=\frac{\alpha}{r^{2}} \frac{d}{d \theta}\left(-\alpha \frac{d u}{d \theta}\right)=-\frac{\alpha^{2}}{r^{2}} \frac{d^{2} u}{d \theta^{2}}\label{eq:ddotr_equation06} \end{eqnarray} By replacing \(\ddot r\) in \eqref{eq:dynamicequ14} with \eqref{eq:ddotr_equation06}, we get: \begin{eqnarray} -\frac{\alpha^{2}}{r^{2}} \frac{d^{2} u}{d \theta^{2}}-\frac{\alpha^{2}}{r^{3}}+\frac{G M}{r^{2}}=0 \quad\text { i.e., }\quad\frac{d^{2} u}{d \theta^{2}}+u=\frac{G M}{\alpha^{2}} \\ \end{eqnarray} This is an easy ordinary differential equation, whose general solution is: \begin{eqnarray} u=A \cos (\theta+\delta)+\frac{G M}{\alpha^{2}}=\frac{1}{r}\label{eq:orbit_quadratic02} \end{eqnarray} where, \(A\) is an arbitrary constant. By introducing the following two parameters, \begin{eqnarray} e:&=&\frac{\alpha^{2}}{G M} A\\ l:&=&\frac{\alpha^{2}}{G M}\\ \end{eqnarray} and we take the origin of angle \(\theta\) so as to make \(\delta=0\) (for example, by assuming \(r\) becomes minimum when \(\theta=0\))), \eqref{eq:orbit_quadratic02} becomes to \begin{eqnarray} r=\frac{l}{1+e \cos \theta}\label{eq:orbit_quadratic04} \end{eqnarray} This is in coincident with a general form of a quadratic curve with \(e\) as the "eccentricity". Now we have arrived at our goal!
Suppose if \(e=0\). Then \eqref{eq:orbit_quadratic04} becomes \(r=l\). It means the orbit is a circle with a radius equal to the constant \(l\).
Suppose if \(0\lt e\lt 1\). Then the denominator of \eqref{eq:orbit_quadratic04} is always a positive value. It is in between \(1-e\) (when \(\theta=\pi\)) and \(1+e\) (when \(\theta=0\)). Therefore, the maximum and the minimum of \(r\) are: \begin{eqnarray} r_\text{max}=\frac{l}{1-e}\quad \text{when } \theta=\pi\\ r_\text{min}=\frac{l}{1+e}\quad \text{when } \theta=0 \end{eqnarray} The satellite moves always in a distance between these two limits. Such curve is (intuitively) an allipse. These two limits are on a common straight line (because \(\theta\) is 0 and \(\pi\)). It means \(r_\text{max}\) + \(r_\text{min}\) is the length of the long axis (major axis). The half of it is called "semi major axis" and described by a constant \(a\) conventionally. Namely,
\begin{eqnarray} a&=&\frac{r_\text{max} + r_\text{min}}{2} = \Bigl(\frac{l}{1-e}+\frac{l}{1+e}\Bigr)/2\nonumber\\ &=&\frac{l}{1-e^2} \end{eqnarray} This is an important formula. Now let's derive the length of the "minor axis", which is perpendicular axis to the major axis. It is the maximum of the \(y\) coordinate of the satellite. If you remember the polar coordinate, you can easily derive \begin{eqnarray} y=r\sin\theta \end{eqnarray} All we need is to find the maximum of this function in terms of \(\theta\). Note that \(r\) depends on \(\theta\). Using \eqref{eq:orbit_quadratic04}, \begin{eqnarray} y=\frac{l\sin\theta}{1+e \cos \theta}\label{eq:y_of_ellipse02} \end{eqnarray} By differentiating it with \(\theta\), we get: \begin{eqnarray} \frac{dy}{d\theta}=...=\frac{l(e+\cos\theta)}{(1+e \cos \theta)^2} \end{eqnarray} It becomes zero when \(e+\cos\theta=0\). In such a circumstance (especially the max of \(y\)), \begin{eqnarray} &&\cos\theta=-e\\ &&\sin\theta=\sqrt{1-\cos^2\theta}=\sqrt{1-e^2} \end{eqnarray} Then \eqref{eq:y_of_ellipse02} takes a valu of \begin{eqnarray} y_{\text {max }}=\frac{l \sqrt{1-e^{2}}}{1-e^{2}}=\frac{l}{\sqrt{1-e^{2}}}=a \sqrt{1-e^{2}} \end{eqnarray} This is the half length of the minor axis, or "semi-minor axis" represented by \(b\) conventionally. Namely, \begin{eqnarray} b=a \sqrt{1-e^{2}} \end{eqnarray} Some peoplemay wonder if \eref{eq:orbit_quadratic04} truely draw an ellipse, even if we have known the length of two main axis. Then I suggest them to try deriving \eref{eq:orbit_quadratic04} by starting from an assumption that the orbit is an ellipse with given values of \(a\) and \(b). Here is a summary of how to do it:
Because the focal point F is at the origin of the coordinate, we need to obtain the position of the center of the ellipse. The distance between the focus and the center is (please consider why by yourself) \((r_\text{max}-r_\text{min})/2=ae\). The ellipse's center is at left with this distance from the origin. Then the equation of the ellipse should be \begin{eqnarray} \frac{(x+ae)^2}{a^2}+\frac{y^2}{b^2}=1 \end{eqnarray} You reform this equation so as to make it only with \(r, e, l, \theta\) by using \begin{eqnarray*} &&x=r\cos\theta\\ &&y=r\sin\theta\\ &&b=a \sqrt{1-e^{2}}\\ &&a=l/(1-e^2) \end{eqnarray*} all of which are already given formula. Then a little complicated calculation (but not very difficult) will derive \eqref{eq:orbit_quadratic04}.
According to the Kepler's second law (or conservation of angular momentum in Newtonian dynamics), the area swept by the line between the focus (the gravity source) and the satellite per unit time (i.e. "aerial speed") is a constant. It is \(r^2\dot\theta/2\). Using \eqref{eq:def_alpha}, it is \(\alpha/2\). (Truely a constant!)
Then the time \(T\) the satellite takes to evolve one time around is, the area of the ellipse \(\pi ab\) divided by this constant aerial speed \(\alpha/2\). Namely, \begin{eqnarray} T=\frac{\pi ab}{\alpha/2} =\frac{\pi a^2\sqrt{1-e^2}}{\sqrt{lGM}/2} =2\pi a^2\sqrt{\frac{1-e^2}{l}}\frac{1}{\sqrt{GM}} =\frac{2\pi a^{3/2}}{\sqrt{GM}} \end{eqnarray} This is a very important formula, which shows the Kepler's third law: The period of the satellite is proportional to (semi major axis)^(3/2). It does not explicitly depend on the minor axis.
But this formula is still not enough. We further want to know at any time where is the satellite on the orbit. We use the following figure. We assume that the satellite is at the point A in the figure at \(t=0\). We represent the area of a shape X as A(X). The time when the satellite arrive at P is determined by A(the fan-like shape FPA), divided by the areal speed.
A(the fan-like shape FPA) is A(the fan-like shape FP'A) \(\times b/a\).
A(the fan-like shape FP'A) is A(fan OP'A) - A(OP'F).
A(fan OP'A) is \(a^2E/2\).
A(OP'F) is \((a\sin E)\times ea /2\).
Therefore, A(the fan-like shape FPA) is \(ab(E-e\sin E)/2\).
By combining all these together, we obtain: \begin{eqnarray} t&=&\frac{ab(E-e\sin E)/2}{\alpha/2} =\frac{\pi ab}{\alpha/2}\frac{E-e\sin E}{2\pi}\nonumber\\ &=&\frac{T}{2\pi}(E-e\sin E)\label{eq:KeplerEq05} \end{eqnarray} Here we define a variable called "mean anomaly" \(M\) as: \begin{eqnarray} M:=\frac{2\pi t}{T}\label{eq:mean_anomaly04} \end{eqnarray} Then \eqref{eq:KeplerEq05} becomes: \begin{eqnarray} M = E-e\sin E\label{eq:KeplerEq08} \end{eqnarray} This is another important formula, called Kepler's equation.
If you want to know the position of the satellite at time \(t\), the following procedures will give you the answer:
