Introduction:

This article on Oribital Mechanics-I is written based on excerpts taken from “Chapter-13: Gravitaion” of my unpublished book on “The Fundamentals of Newtonian Mechanics”. This is basically a kind of review paper addressing the fundamental mathematical tools that are required to understand the orbital mechanics. It is intended to inspire my present and past students and professor colleagues who are interested to develop the next generation “Deep Space Communication and Exploration of Solar System through Inter-Lagrangian Data Relay Satellite Constellation”, which was presentated by my students at iCubeSat 2019, the 8th Interplanetary CubeSat Workshop at the Politecnico di Milano, Milan, italy on 28-29 May 2019.

Rocket and space technology is based on the application of orbital mechanics or flight mechanics. Orbital mechanics deals with the motions of artificial satellites and space vehicles under the influence of gravity, atmospheric drag, and thrust. Orbital mechanics is a modern branch of celestial mechanics which is the study of the motions of natural celestial bodies such as the star, galaxies, moon and planets. Isaac Newton (1642-1727), a 17th century mathematician who put forward the laws of motion and formulated the law of universal gravitation is the theoretical founding principle of orbital mechanics. The engineering applications of orbital mechanics include ascent trajectories, reentry and landing, numerical modeling and simulations, lunar, interplanetary trajectories, communication sattelites, space stations,space shuttles and inter-planetary communication.

Newtonian Gravity

Newtonian mechanics was initially developed in order to account for the motion of the Planets around the Sun. All matter exerts a force, called gravity, that pulls all other matter towards its center. The strength of the force depends on the mass of the object. The Sun is about 328 times heavier than the Earth and as a result it has more gravity than Earth. As an example Jupiter is 318 times heavier than Earth, and its gravity is \(22.9\space m/s^2\) whereas Earth’s gravity on its surface is \(9.8~m/s^2\). Objects far from the Sun won’t be influenced by its gravity. The force that keeps the Planets in their respective orbits around the Sun is called gravitational force and was first described by Isaac Newton in his Philosophiae Naturalis Principia Mathematica (“Mathematical Principles of Natural Philosophy,”) published in 1687. According to Newton, any two point mass objects (or spherically symmetric objects of finite extent) exert a force of attraction on one another. This force points along the line of centers joining the objects, is directly proportional to the product of their masses, and inversely proportional to the square of the distance between them.

Let’s assume that the first object is the Sun, with mass \(M\), and is located at the origin of our coordinate system. Let’s assume that the second object be a planet, of mass \(m\), which is located at position vector \({\vec r}\). The gravitational force exerted on the planet by the Sun is thus written

\[\begin{equation} \tag{13-1} \vec F = - \frac{G\,M\,m}{r^3}\,{\vec r}. \end{equation}\]

The constant of proportionality, \(G\), is called the universal gravitational constant, and its value in SI unit is given by:

\[\begin{equation} \tag{13-2} G = 6.67300\times 10^{-11}\space \frac{m^3}{kg-s^2}. \end{equation}\]

We shall assume that the Sun is so much more massive than the planet in question that the force exerted by the planet does not cause the Sun’s position to shift appreciably. Hence, the Sun will always remain at the origin of our coordinate system. Likewise, we shall neglect the gravitational forces exerted on our planet by the other planets in the Solar System, compared to the much larger gravitational force exerted by the Sun.

According to Eq. 1, the gravitational force acting on an object is directly proportional to its inertial mass. This question perplexed physicists for many years, and was only answered when Albert Einstein published his general theory of relativity in 1916. According to Einstein, inertial mass acts as a sort of gravitational charge since it is impossible to distinguish an acceleration produced by a gravitational field from an apparent acceleration generated by observing in a non-inertial reference frame. The assumption that these two types of acceleration are indistinguishable leads directly to all of the strange predictions of general relativity such as clocks in different gravitational potentials run at different rates, time-space curvature etc.

Applying Newton’s second law in Eq. 1, the equation of motion of our planet takes the following form

\[\begin{equation} \tag{13-3} m\frac{d^2{\vec r}}{dt^2} = -\frac{GmM}{r^3}\,{\vec r}. \end{equation}\]

\[\begin{equation} \tag{13-4} \frac{d^2{\vec r}}{dt^2} = -\frac{G\,M}{r^3}\,{\vec r}. \end{equation}\]

Note that the planetary mass, \(m\), has canceled out from both sides of the above equation.

Gravitation Near Earth’s Surface:

Assuming Earth is a uniform sphere of mass M. Please don’t get confused with capital M, being the mass of earth now. In the last paragraph we considered M being the mass of Sum and m was the mass of a planet. The magnitude of the gravitational force from Earth on a particle of mass m, located outside Earth at a distance r from Earth’s center, is then given by

\[\begin{equation} \tag{13-5} {F} = \frac{G\,M\,m}{r^2}\ \end{equation}\]

all objects fall toward the center of Earth due to gravitational field around the Earth . The result of gravitational force \(\vec F\) and its associated acceleration we call the gravitational acceleration \(\vec a_g\). Newton’s second law tells us that magnitudes F and \(a_g\) are related by

\[\begin{equation} \tag{13-6} {F} = ma_g \end{equation}\]

\[\begin{equation} \tag{13-7} a_g = \frac{G\,M}{r^2}\ \end{equation}\]

Table 13-1 shows values of \(a_g\) computed for various altitudes above Earth’s surface. Notice that \(a_g\) is significant even at 400 km.

Table 13-1 Variation of \(a_g\) with Altitude

Altitude (km) \(a_g (m/s^2)\) Altitude Example
0 9.83 Mean Earth surface
8.8 9.80 Mt. Everest
36.6 9.71 Highest crewed balloon
400 8.70 Space shutte orbit
35700 0.225 Communications satellite

In Chapter 5, of my book on “The Fundamentals of Newtonian Mechanics” we made an assumption that Earth is a non-accelerating inertial frame by neglecting its rotation on its own axis. This assumption led us to assume that the free-fall acceleration \(g\) of a particle is the same as the particle’s gravitational acceleration \(a_g\). Also, we assumed that \(g\) has the constant value \(9.8 ~\frac {m}{s_2}\) any place on Earth’s surface. However, any \(g\) value measured at a given location will differ from the \(a_g\) value calculated with Eq. 13-7 for that location for three reasons: (1) Earth’s mass is not distributed uniformly (Fig. 13-1), (2) Earth is not a perfect sphere, and (3) Earth rotates about its own axis through the poles.

Moreover, because \(g\) differs from \(a_g\), as a result the measured weight \(mg\) of a particle differs from the magnitude of the gravitational force on the particle as given by Eq. 5. Let us now examine those reasons.

1. Non-uniform Distribution of Earth’s Mass. The density (mass per unit volume) of Earth varies radially as shown in Fig. 1, as a result g varies from region to region over the surface.

2. Earth is not a sphere. Earth is slightly ellipsoidal in shape. Its equatorial radius (from its center point out to the equator) is greater by 21 km than its polar radius (from its center point out to either north or south pole). On in other words a point at the poles is closer to the dense core of Earth than is a point on the equator. The free-fall acceleration g slightly more than the equator (see Eq. 13-7).

3. Earth is rotating. Earth is rotating with a period of T = 24 hrs about an axis of rotation passing through the poles. A rotating Earth has maximum centripetal acceleration or a centripetal net force directed toward the Earth’s center at the equator and minimum at the poles.

The effect of rotation on g can be easily calculated by considering a free-body diagram for the crate as shown in Figure 13-2. As can be seen in Figure 13-2, two forces acting on the crate, one is the normal force \(\vec F_N\) directed radially upward and the other one is The gravitational force \(m\vec g\) directed radially towards the center of the Earth.

Using Newton’s second law for forces along the radial direction

\[\begin{equation} \tag{13-8} F_{net,r} = ma_r \end{equation}\]

\[\begin{equation} \tag{13-9} F_N -ma_g = m(-\omega^2R) \end{equation}\]

The magnitude \(F_N\) of the normal force is equal to the weight mg read on the scale. Substituting \(F_N = mg\) in Eq. 23-9 gives us

\[\begin{equation} \tag{13-10} mg = ma_g -m(\omega^2R) \end{equation}\]

\[\begin{equation} \tag{13-11} g = a_g -(\omega^2R) \end{equation}\]

where \(\omega=\frac{2\pi}{T}\), with \(T=24~h\) is the period of rotation of Earth on its own axis. The measured free-fall acceleration is less than the gravitational acceleration because the of Earth’s rotation.

Considering \(T=24~h\) and radius of Earth \(R=6.37\times 10^6\space m\), he value of \(\omega^2R=\frac{4pi^2R}{T^2}=0.03369~m/s^2\), and which is small compared to \(g=9.8~m/s^2.\) This little centripetal acceleration at the equator due to Earth’s rotation on its own axis make the Earth slightly non-inertial.

Gravitational Potential Energy:

We want to find the gravitational potential energy U of the base ball at point P along its path, at radial distance R from Earth’s center as shown in Fig. 13-3.

To find the gravitational potential energy U, we first find the work W done on the ball by the gravitational force as the ball travels from point P to a great (infinite) distance from Earth. Because the gravitational force (r) is a variable force (its magnitude depends on r). In vector notation, we can write

\[\begin{equation} \tag{13-12} W = \int_R^{\infty}\vec F(r).d\vec r = \int_R^{\infty}F(r)dr\cos\phi = \int_R^{\infty}F(r)dr\cos180^0=-\int_R^{\infty}F(r)dr \end{equation}\]

substituting Eq. 13-5 in Eq. 13-12, we get the final expression for W, as follows:

\[\begin{equation} \tag{13-13} W = -\int_R^{\infty}F(r)dr=-{GMm}\int_R^{\infty}\frac{1}{r^2}dr=-\frac{GMm}{R} \end{equation}\]

where M is the mass of Earth and m is the mass of the baseball taken into consideration in this example. We know from Eq. 8-1 (\(\Delta U = -W\)) tells us that we can also write that work in terms of potential energies as

\[U_{\infty}-U = -W\] \[U = W=-\frac{GMm}{R},\] where we assumed the gravitational potential engery, \(U_{\infty}=0\), at infinity.

Escape Speed:

If a projectile is thrown upward with a certain minimum initial speed that will cause it to move upward forever, theoretically coming to rest only at infinity. This minimum initial speed is called the (Earth) escape speed.

Applying the principle of conservation of energy, the total energy at the planet’s surface must be zero

\[K+U = \frac{1}{2}mv^2+\left(-\frac{GMm}{R}\right)=0\]

\[\begin{equation} \tag{13-14} v=\sqrt\frac{2GM}{R} \end{equation}\]

Note that v does not depend on the direction in which a projectile is fired from a planet. However, attaining that speed is easier if the projectile is fired in the direction the launch site is moving as the planet rotates about its axis. For example, rockets are launched eastward at Cape Canaveral to take advantage of the Cape’s eastward speed of 1500 km/h due to Earth’s rotation.

Eq. 13-14 can be used to find the escape speeds for various astronomical objects, taking into consideration of the masses M and radii R of the respective objects under consideration. Table 13-2 shows some escape speeds.

Table 13-2 Some Escape Speeds

Body Mass (kg) Radius (m) Escape speed (km/s)
\(Ceres^a\) \(1.17\times 10^{21}\) \(3.8\times 10^5\) 0.64
\(Earth's moon^a\) \(7.36\times 10^{22}\) \(1.74\times 10^{6}\) 2.38
Earth \(5.98\times 10^{24}\) \(6.37\times 10^{6}\) 11.2
Jupiter \(1.90\times 10^{27}\) \(7.15\times 10^{7}\) 59.5
Sun \(1.99\times 10^{30}\) \(6.96\times 10^{8}\) 618
\(Sirius B^b\) \(2\times 10^{30}\) \(1\times 10^{7}\) 5200
\(Neutron star^c\) \(2\times 10^{30}\) \(1\times 10^{4}\) \(2\times 10^5\)

\(^a\) The most massive of the asteroids.

\(^b\) A white dwarf (a star in a final stage of evolution) that is a companion of the bright star Sirius.

\(^c\) The collapsed core of a star that remains after that star has exploded in a supernova event.

Conservative Force:

We have seen the force of gravity is a conservative force in Eq. 5-21. Hence, we can rewrite the gravitational force (13-1) in terms of Eq.5-21

\[\begin{equation} \tag{13-15} {\vec F} = - \vec \nabla U(\vec r), \end{equation}\]

where the potential energy, \(U({\vec r})\), of our planet in the Sun’s gravitational field takes the form

\[\begin{equation} \tag{13-16} U({\vec r}) = - \frac{G\,M\,m}{r^2}{\vec r}. \end{equation}\]

It follows that the total energy of our planet is a conserved or constant quantity. In other words,

\[\begin{equation} \tag{13-17} {\ E} = {\ K} + {\ U} = \frac{mv^2}{2} - \frac{G\,Mm}{r} \end{equation}\]

or in other words, total energy is the sum of kinetic energy (K) and potential energy (U) is constant in time. Here, \({E}\) is the planet’s total energy, and \({\vec v} = d{\vec r}/dt\) is the linear velocity.

Gravity is also a central force and the angular momentum of our planet is a conserved quantity. In other words,

\[\begin{equation} \tag{13-18} {\vec l} = {\vec r}\times {\vec p}, \end{equation}\]

where \(\vec p = m\vec v\) is the linear momentum, \(\vec l\) is the planet’s angular momentum, is constant in time. Taking the scalar product of the above equation with \({\vec r}\), we obtain

\[\begin{equation} \tag{13-19} {\vec l}\cdot{\vec r} = 0. \end{equation}\]

This is the equation of a plane which passes through the origin, and whose normal is parallel to \({\vec l}\). Since \({\vec l}\) is a constant vector, it always points in the same direction. We, therefore, conclude that the motion of our planet is two-dimensional in nature: i.e., it is confined to some fixed plane which passes through the origin. Without loss of generality, we can let this plane coincide with the \(x\)-\(y\) plane.

Polar Coordinate System:

Polar coordinate system is shown in Figure 13-4 with respect to \(x\)-\(y\) cartesian coordinate system. We can determine the instantaneous position of our planet in the \(x\)-\(y\) plane in terms of standard Cartesian coordinates, (\(x\), \(y\)), or polar coordinates, (\(r\), \(\theta\)), as can be seen in Figure 13-4.

The magnitude of the position vector \(\vec r\) and direction are, \(r=\sqrt{x^2+y^2}\) and \(\theta=\tan^{-1}(y/x)\) respectively. The instantaneous position of the planet can be defined in terms of two unit vectors, \({\hat e}_r\equiv {\hat r}/r\) and \({\hat e}_\theta\equiv {\hat e}_z\times {\hat e}_r\). The position unit vector points radially outward from the origin, whereas the second is normal to the first, in the direction of increasing \(\theta\). The Cartesian components of \({\hat e}_r\) and \({\hat e}_\theta\) are

\[\begin{equation} \tag{13-20} {\hat e}_r = (\cos\theta,\sin\theta), \end{equation}\]

\[\begin{equation} \tag{13-21} {\hat e}_\theta = (-\sin\theta,\cos\theta). \end{equation}\]

The position vector \(\vec r\) of the planet in terms of polar coordinate is shown in Figure 13-5. In Figure 13-5, the semi-major axis, perihelion (closest to the Sun) and aphelion (farest from the Sun) distances are represented by a, \(R_p\) and \(R_a\) respectively.

\[\begin{equation} \tag{13-22} \vec r = r\hat e_r \end{equation}\]

Thus, the planet’s linear velocity \(\vec v\) becomes

\[\begin{equation} \tag{13-23} {\vec v} = \frac{d{\vec r}}{dt} = \dot{r}\,{\hat e}_r + r\,\dot{\hat e}_r \end{equation}\]

where \(\dot{~}\) is shorthand for \(d/dt\). Note that \({\hat e}_r\) has a non-zero time-derivative (unlike a Cartesian unit vector) because its direction changes as the planet moves around. As is easily demonstrated, from differentiating Eq. 13-10 with respect to time,

\[\begin{equation} \tag{13-24} \dot{\hat e}_r = \dot{\theta}\,(-\sin\theta,\,\cos\theta) = \dot{\theta}\,\,{\hat e}_\theta. \end{equation}\]

Thus, \[\begin{equation} \tag{13-25} {\vec v} = \dot{r}\,\,{\hat e}_r + r\,\dot{\theta}\,\,{\hat e}_\theta. \end{equation}\]

Now, the planet’s acceleration is written

\[\begin{equation} \tag{13-26} {\vec a} = \frac{d{\vec v}}{dt} = \frac{d^2{\vec r}}{dt^2} = \ddot r \hat e_r + \dot r \dot{\hat e}_r +(\dot r\dot\theta+r\ddot\theta){\hat e}_\theta + r\,\dot{\theta}\,\,\dot{\hat e}_\theta. \end{equation}\]

Again, \({\hat e}_\theta\) has a non-zero time-derivative because its direction changes as the planet moves around. Differentiation of Eq. 13-11 with respect to time yields

\[\begin{equation} \tag{13-27} \dot{\hat e}_\theta = \dot{\theta}\,(-\cos\theta,\,-\sin\theta) = - \dot{\theta}\,{\hat e}_r. \end{equation}\]

Hence,

\[\begin{equation} \tag{13-28} {\vec a} = (\ddot{r}-r\,\dot{\theta}^{\,2})\,{\hat e}_r + (r\,\ddot{\theta} + 2\,\dot{r}\,\dot{\theta})\,{\hat e}_\theta. \end{equation}\]

It follows that the equation of motion of our planet, (Eq. 13-4), can be written

\[\begin{equation} \tag{13-29} {\vec a} = (\ddot{r}-r\,\dot{\theta}^{\,2})\,{\hat e}_r + (r\,\ddot{\theta} + 2\,\dot{r}\,\dot{\theta})\,{\hat e}_\theta = - \frac{G\,M}{r^2}\,{\hat e}_r. \end{equation}\]

Since \({\hat e}_r\) and \({\hat e}_\theta\) are mutually orthogonal, we can separately equate the coefficients of both, in the above equation, to give a radial equation of motion,

\[\begin{equation} \tag{13-30} \ddot{r}-r\,\dot{\theta}^{\,2} = - \frac{G\,M}{r^2}, \end{equation}\]

and a tangential equation of motion,

\[\begin{equation} \tag{13-31} r\,\ddot{\theta} + 2\,\dot{r}\,\dot{\theta} = 0. \end{equation}\]

Kepler’s Laws

Planetary motion in the Solar System was first correctly described by astronomer Johannes Kepler in his published works between 1609 and 1619. He described the motion of the Planets in the following three simple laws:

Now let’s see how Kepler’s laws can be explained and derived from Newton’s laws of motion.

Kepler’s First Law:

All planetary orbits are elliptical which are confocal with the Sun (i.e., the Sun lies at one of the focii of each ellipse).

Multiplying both sides of Eq.13-31 by r we get,

\[\begin{equation} \tag{13-32} r^2\,\ddot{\theta} + 2r\,\dot{r}\,\dot{\theta} = 0. \end{equation}\]

Eq. 13-32, can also be expressed as

\[\begin{equation} \tag{13-33} \frac{d{(r^2\dot\theta)}}{dt}=0 \end{equation}\]

Eq. 13-33, implies that

\[\begin{equation} \tag{13-34} l = r^2\dot\theta = constant \end{equation}\]

We already proved this in Eq. 5-39, using Newton’s Third Law of Motion. We have discovered that \(l\) is constant in time and is the angular momentum of our planet. In discussing Newtoian mechanics in Chapter-5 we already mentioned that this is the case because gravity is a central force.

Our planet’s radial equation of motion, (Eq. 13-30), can be combined with Eq.13-34 to give

\[\begin{equation} \tag{13-35} \ddot{r} -\frac{l^2}{r^3}= - \frac{G\,M}{r^2} \end{equation}\]

Suppose that \(r = u^{-1}\). It follows that

\[\begin{equation} \tag{13-36} \dot{r} = - \frac{\dot{u}}{u^2} = - r^2\,\frac{du}{d\theta}\,\frac{d\theta}{dt} = - l\,\frac{du}{d\theta} \end{equation}\]

Similarly,

\[\begin{equation} \tag{13-37} \ddot{r} = - l \,\frac{d^2 u}{d\theta^2}\,\dot{\theta} = - u^2\,l^2\,\frac{d^2 u}{d\theta^2} \end{equation}\]

Hence, Eq. 13-35 can be written in the linear form

\[\begin{equation} \tag{13-38} \frac{d^2 u}{d\theta^2} + u = \frac{G\,M}{l^2} \end{equation}\]

The general solution to the above equation is

\[\begin{equation} \tag{13-39} u(\theta) = \frac{G\,M}{l^2}\left[1 - e\,\cos(\theta-\theta_0)\right] \end{equation}\]

where \(e\) and \(\theta_0\) are arbitrary constants. Without loss of generality, we can set \(\theta_0=0\) by rotating our coordinate system about the \(z\)-axis. Thus, we obtain

\[\begin{equation} \tag{13-40} r(\theta) = \frac{r_c}{1 - e\,\cos\theta} \end{equation}\]

where

\[\begin{equation} \tag{13-41} r_c = \frac{l^2}{G\,M} \end{equation}\]

We find Eq. 13-39 as the equation of a conic section (see Eq. 13-35) which is confocal with the origin (i.e., with the Sun). Specifically, for \(e<1\), Eq. 13-39 is the equation of an ellipse which is confocal with the Sun. Thus, the orbit of our planet around the Sun in a confocal ellipse–this is Kepler’s first law of planetary motion. Of course, a planet cannot have a parabolic or a hyperbolic orbit, since such orbits are only appropriate to objects which are ultimately able to escape from the Sun’s gravitational field.

Conic Sections:

The ellipse, the parabola, and the hyperbola are collectively known as conic sections, since these three types of curve can be obtained by taking various different plane sections of a right cone. It turns out that the possible solutions of Eq. 13-30 and 13-31 are all conic sections.

Parabola:

A parabola which is aligned along the \(+x\)-axis, and passes through the origin (see Fig. 13-6).

\[\begin{equation} \tag{13-42} y^2 - 4ax = 0 \end{equation}\]

where \(a>0\).

In the case of the parabola the eccentricity \(e = 1\) For ellipses, the eccentricity \(0 < e < 1\) and it is defined as \[e = \sqrt{\left(1-\frac{b^2}{a^2}\right)}\]

The foci are \(\pm ae,0\) and the directrix are \(x = \pm \frac{a}{e}\)$$

Ellipse:

An ellipse, centered on the origin, of major radius \(a\) and minor radius \(b\), which are aligned along the \(x\)- and \(y\)-axes, respectively (see Fig. 13-2), satisfies the following well-known equation:

\[\begin{equation} \tag{13-43} \frac{x^2}{a^2} + \frac{y^2}{b^2} = 1,\space where\space a>b>0. \end{equation}\]

In case of circles a = b. A circle in standard position with an ellipse is shown in Figure 13-7.

Hyperbola:

A hyperbola which is aligned along the \(+x\)-axis, and whose asymptotes intersect at the origin (see Figure 13-8).

\[\begin{equation} \tag{13-44} \frac{x^2}{a^2} - \frac{y^2}{b^2} = 1 \end{equation}\]

The eccentricity of hyperbola, \[e = \sqrt{\left(1-\frac{b^2}{a^2}\right)}\] The foci and directrices are the same as the ellipsis and the asymptotes can be found at \(y=\pm \frac{b}{a}x\).

In Figure-13-8, \(a\) is the distance of closest approach to the origin. The asymptotes subtend an angle \(\phi=\tan^{-1}(b/a)\) with the \(x\)-axis. All conic sections can be described in terms of eccentricity. Satellite orbits are any of the four conic sections. All conic sections are have semi-major axis and the energy. The Table 13-3 shows the type of conic sections and their relationships with eccentricity, semi-major axis, and energy.

Table 13-3 Conic Section Parameters

Conic Section Eccentricity, e Semi-major Axis Energy
Circle 0 a=Radius E<0
Ellipse 0 < e < 1 a>0 E<0
Parabola e = 1 a=Infinity E=0
Hyperbola e > 1 a<0 E>0

In summarizing the discussion on conic section, or just conic, is a curve formed by passing a plane through a right circular cone. As shown in Figure 13-9, the angular orientation of the plane relative to the cone determines whether the conic section is a circle, ellipse, parabola, or hyperbola. The circle and the ellipse arise when the intersection of cone and plane is a bounded curve. The circle is a special case of the ellipse in which the plane is perpendicular to the axis of the cone. If the plane is parallel to a generator line of the cone, the conic is called a parabola. If the intersection is an unbounded curve and the plane is not parallel to a generator line of the cone, the figure is a hyperbola. In the latter case the plane will intersect both halves of the cone, producing two separate curves.