Orbital Trajectory

To derive the orbit equation we must first begin by describing the energy of the reduced mass system orbiting the system's center of mass. We know that the kinetic energy and the gravitational potential of the system, so we can express the total energy:

$$E=\frac{1}{2}\mu v^2-\frac{Gm_1m_2}{r},$$

where \mu is the reduced mass in \(\text{kg}\) and \(v\) is the relative velocity between the two bodies.

The velocity can be defined in polar coordinates as,

\begin{align*}\vec{v}&=v_\text{r} \hat{r} +v_\theta \hat{\theta},\\v&=|\vec{v}|,\\v&=|\frac{\text{d}\vec{r}}{\text{d}t}|,\end{align*}

where \(v_\text{r}=\frac{\text{d}r}{\text{d}t}\) and \(v_\theta=r\left(\frac{\text{d}\theta}{\text{d}t}\right)\).

The energy of the system is now:

$$E=\frac{1}{2}\mu\left[\left(\frac{\text{d}r}{\text{d}t}\right)^2+\left(r\frac{\text{d}\theta}{\text{d}t}\right)^2\right]-\frac{Gm_1m_2}{r}.$$

The angular momentum is given by

$$\begin{align*}\vec{l}&=\vec{r}\times\mu\vec{v},\\\vec{l}&=r \hat{r}\times\mu\left(v_\text{r} \hat{r}+v_\theta \hat{\theta}\right),\\\vec{l}&=r\mu v_\theta \hat{k},\\\vec{l}&=\mu r^2\frac{\text{d}\theta}{\text{d}t} \hat{k},\\l&=\mu r^2\frac{\text{d}\theta}{\text{d}t},\\\frac{\text{d}\theta}{\text{d}t}&=\frac{l}{\mu r^2}.\end{align*}$$

Now we rewrite the total energy of the system,

$$E=\frac{1}{2}\mu\left(\frac{\text{d}r}{\text{d}t}\right)^2+\frac{1}{2}\frac{l^2}{\mu r^2}-\frac{Gm_1m_2}{r}.$$

We can write it in terms of \(\frac{\text{d}r}{\text{d}t}\):

$$\frac{\text{d}r}{\text{d}t}=\sqrt{\frac{2}{\mu}}\left(E-\frac{1}{2}\frac{l^2}{\mu r^2}+\frac{Gm_1m_2}{r}\right)^{\frac{1}{2}}$$

Now we divide \(\frac{\text{d}\theta}{\text{d}t}\) by \(\frac{\text{d}r}{\text{d}t}\) to obtain \(\frac{\text{d}\theta}{\text{d}r}\), an expression that relates the distance \(r\) to the angle \(\theta\). It gives the orbit equation in differential form:

\begin{align*}\frac{\text{d}\theta}{\text{d}r}&=\frac{\frac{\text{d}\theta}{\text{d}t}}{\frac{\text{d}r}{\text{d}t}},\\\frac{\text{d}\theta}{\text{d}r}&=\frac{l}{\sqrt{2\mu}}\frac{\left(\frac{1}{r^2}\right)}{\left(E-\frac{l^2}{2\mu r^2}+\frac{Gm_1m_2}{r}\right)^{\frac{1}{2}}},\\\text{d}\theta&=\frac{l}{\sqrt{2\mu}}\frac{\left(\frac{1}{r^2}\right)}{\left(E-\frac{l^2}{2\mu r^2}+\frac{Gm_1m_2}{r}\right)^{\frac{1}{2}}}\text{d}r.\end{align*}

Before we do any integration, we need to do some substitutions. We do the substitution \(u=\frac{1}{r}\), \(\text{d}u=-\left(\frac{1}{r^2}\right)\text{d}r\,\):

$$\text{d}\theta=-\frac{l}{\sqrt{2\mu}}\frac{\text{d}u}{\left(E-\frac{l^2}{2\mu}u^2+Gm_1m_2u\right)^{\frac{1}{2}}}.$$

Next, we multiply and divide the right hand side of the equation by \(\frac{\sqrt{2\mu}}{l}\):

$$\begin{align*}\text{d}\theta&=-\frac{\text{d}u}{\left(\frac{2\mu E}{l^2}-u^2+2\left(\frac{\mu Gm_1m_2}{l^2}\right)u\right)^{\frac{1}{2}}},\\\text{d}\theta&=-\frac{\text{d}u}{\left(\frac{2\mu E}{l^2}-u^2+\frac{2u}{r_0}\right)^{\frac{1}{2}}},\end{align*}$$

where we have defined \(r_0=\frac{l^2}{\mu Gm_1m_2}\).

Now we add and subtract \(\frac{1}{{r_0}^2}\) inside the parenthesis with the square root:

\begin{align*}\text{d}\theta&=-\frac{\text{d}u}{\left(\frac{2\mu E}{l^2}+\frac{1}{{r_0}^2}-u^2+\frac{2u}{r_0}-\frac{1}{{r_0}^2}\right)^{\frac{1}{2}}},\\\text{d}\theta&=-\frac{\text{d}u}{\left(\frac{2\mu E}{l^2}+\frac{1}{{r_0}^2}-\left(u-\frac{1}{r_0}\right)^2\right)^{\frac{1}{2}}},\\\text{d}\theta&=-\frac{r_0\text{d}u}{\left(\frac{2\mu E{r_0}^2}{l^2}+1-\left(r_0 u-1\right)^{2}\right)^{\frac{1}{2}}},\end{align*}

where we define the eccentricity as \(\epsilon=\sqrt{\frac{2\mu E{r_0}^2}{l^2}+1}\). This dimensionless quantity is responsible for the shape of the orbit. This is better discussed in the article Orbital Trajectories.

We rewrite our equation in terms of the eccentricity:

$$\text{d}\theta=-\frac{r_0\text{d}u}{\left({\epsilon}^2-\left(r_0 u-1\right)^{2}\right)^{\frac{1}{2}}}.$$

Finally, we do the last substitution before solving the integral, \(r_0u-1=\epsilon \cos{\alpha}\) and \(r_0 \text{d}u=-\epsilon\sin{\alpha}\text{d}\alpha\):

\begin{align*}\text{d}\theta&=-\frac{-\epsilon\sin{\alpha}\text{d}\alpha}{\left({\epsilon}^2-{\epsilon}^2{\cos^2{\alpha}}\right)^{\frac{1}{2}}},\\\text{d}\theta&=-\frac{-\bcancel{\epsilon}\sin{\alpha}\text{d}\alpha}{\bcancel{\epsilon}\left(1-\cos^2{\alpha}\right)^{\frac{1}{2}}},\\\text{d}\theta&=\frac{\sin{\alpha}\text{d}\alpha}{\left(sin^2{\alpha}\right)^{\frac{1}{2}}},\\\text{d}\theta&=\frac{\bcancel{\sin{\alpha}}\text{d}\alpha}{\bcancel{sin{\alpha}}}\,\\\theta&=\int{\text{d}\alpha},\\\theta&=\alpha + \text{constant}.\end{align*}

We already know that \(r=\frac{1}{u}\). If we choose the constant to be zero, we find that:

\begin{align*}r_0u-1&=\epsilon\cos{\alpha},\\r_0u-1&=\epsilon\cos{\theta},\\u&=\frac{1-\epsilon\cos{\theta}}{r_0}.\end{align*}

Our final expression for the orbit equation is:

\begin{align*}r&=\frac{1}{u},\\r&=\frac{r_0}{1+\epsilon\cos{\theta}}.\end{align*}

Alternately, if we choose the constant to be \(\pi\) we get the same orbit equation but with the vertical axis reflected,

$$r=\frac{r_0}{1-\epsilon\cos{\theta}}.$$

We need two coordinates to track the position of an orbiting satellite in the sky, the azimuth and elevation angles. The azimuth angle is an angle of \(360^\circ\) degrees of rotation horizontally, while the elevation angle is an angle of \(90^\circ\) degrees of rotation vertically. These angles are measured from the antenna on Earth. We present a diagram below to help you visualize these angles.

Fig. 5 - Diagram to visualize the azimuth and elevation angles and be able to track the position of a satellite in orbit.

The equation to calculate the azimuth and elevation angles varies with the type of orbit. Here we will discuss the equation for satellites in a geostationary orbit. If we want the more general equation, it gets a little bit more complicated as we have to add some factors like the satellite's height above sea level and a time-dependent equation, as not every orbit has the same period as the Earth.

We can express the equation for the azimuth and elevation angles of geostationary satellites as follows:

$$A=180^\circ+\;\tan^{-1}\left(\frac{\tan\left(S-N\right)}{\tan\left(L\right)}\right)$$

$$E=\tan^{-1}\left(\frac{\cos\left(S-N\right)\cos\left(L\right)-0.1512}{\sqrt{1-\cos^2\left(S-N\right)\cos^2\left(L\right)}}\right),$$

where \(S\) is the satellite's longitude in degrees, \({}^\circ\), \(N\) is the longitude of the antenna site in degrees, \({}^\circ\), and \(L\) is the latitude of the antenna site in degrees, \({}^\circ\).