MathJAX

Thursday, 6 July 2017

World Population Distribution by Latitude & Longitude

After seeing the excellent stuff on this topic from Radical Cartography and Datagraver, I wanted to have my own version.

I started with the source data set from SEDAC/CIESIN (requires registration) which is used by both the sites named above. I used the 2015 data set, which the Datagraver page also uses. The data comes in a GeoTIFF file which contains the population count 'pixel' for every latitude/longitude point for every 30 arc seconds (½ arc minutes or 1/120 degree) of latitude & longitude. Since the Earth's circumference is just over 40,000 km, this means that at a latitude λ this 'pixel' is a rectangle (a trapezium if we want to be pedantic, but never mind) that is 926 meters long and 926 cos λ meters wide: at the equator it is a 926m × 926m square. Actually, the data only covers latitudes from 85°N to 60°S which leaves out Antarctica, and it also marks the interior of Greenland as 'no data', but we can live with that!

This turned out to be a big exercise. For starters, the GeoTIFF file is over 400 MB in size. I used QGIS to extract the data from the TIFF file into a 3-column 'longitude, latitude, value' text file:

  1. Start QGIS
  2. 'Project' ▶ 'New'
  3. 'Layer' ▶ 'Add Layer' ▶ 'Add Raster Layer...'
  4. Select the GeoTIFF file and 'Open'
  5. Wait for file to load (may take a while)
  6. 'Raster' ▶ 'Conversion' ▶ 'Translate (Convert Format)...'
  7. Select output file type = 'ASCII Gridded XYZ'
  8. Specify output file, then 'OK'
  9. Wait ... a long time

This resulted in a file which was over 40 GB in size: not a surprise considering that it had (80 + 65) × 120 × 360 × 120 = 751,680,000 data points in it. I found that 'no data' is indicated by the special value -407649103380480, and running a 'fgrep -v' to remove those points gave me a file that was still over 11 GB (with 212,565,537 data points.)

Finally, after running this through a jury-rigged C-program to get the data down to a granularity of 1° instead of 30", I arrived at a 650 kB text file. This I imported into Excel, and got these graphs.

The graph drawn in red is the population distribution by longitude. The one in blue is the distribution by latitude, while the one in green is the distribution by 'folded latitude' which is the absolute value of the latitude or simply the distance from equator.

For anyone who needs the raw data (in 1° intervals), here is the spreadsheet. The raw data is in D3:MZ183, the rest are subtotals etc.

Some statistics:

  • The mean latitude is 22.7°N, while the median latitude—the 'population equator'—is 25.2°N. The highest mode is at 26°N, followed by 23°N.
  • Going by the distance from equator, the mean is at 26.2°. The median is at 25.8°, which means half the population lives within this distance of the equator. The highest mode is at 23°, followed by 26°.
  • 80% of the population lives within 36.9° of the equator, 90% within 43.5°, 95% within 49.6° and 99% within 55.6°. Here is a graph of the cumulative population distribution by distance from the equator:
  • The longitude with the highest population is 77°E, on which Delhi & Bangalore lie. This is followed by 121°E, on which Shanghai & Manila lie.
  • The southern hemisphere is home to only about 13% of the world's population.

Here is my attempt at plotting the population data in a logarithmic scale 'heat-map' for the 1° grid. I've used violet for less than 100, light blue for 100 to 1,000, darker blue for 1,000 to 10,000, green for 10,000-100,000, yellow for 100,000-1,000,000, orange for 1,000,000-10,000,000 and red for more than 10,000,000.


Monday, 15 May 2017

A little bit of satellite geometry

This is a look at the geometry and trigonometry of satellites and their orbits. To a large extent, this post falls back to describing circular orbits around a spherical Earth, but there are exceptions that (I hope) makes things interesting.

Looking up at the sky

We have an observer at latitude $\varphi$ and longitude $\lambda$. Also, we have a satellite in the sky at an altitude $h$ directly above a 'subsatellite point' which is located at latitude $\varphi_s$ and longitude $\lambda_s = \lambda + \Delta\lambda$. Which way must the observer look to see that satellite?

By definition, latitudes $\varphi$ and $\varphi_s$ are between $-90^\circ$ and $+90^\circ$ (N: +ve, S: -ve) and the $\Delta\lambda$ is between $-180^\circ$ and $+180^\circ$ (E: +ve, W: -ve).

The direction towards which the observer must look can be expressed in terms of two angles.

  1. The azimuth $\psi$ (also known as bearing) which is an angle from $-180^\circ$ clockwise to $+180^\circ$ with $0^\circ$ being due North.
  2. The elevation $\beta$ (also known as the look angle when we are talking about satellites) which is an angle from $-90^\circ$ (nadir) upwards to $+90^\circ$ (zenith) with $0^\circ$ being towards the horizon.

For negative values of $\beta$ the satellite is not visible because it is below the horizon, so for our purposes we can take $\beta$ to be between $0^\circ$ and $+90^\circ$.

The diagram above illustrates what we mean by $\psi$ and $\beta$. It also introduces the distance $d$ between the observer and the satellite, and the distance $D$ along the ground between the observer and the subsatellite point. Of course, unlike this simplified illustration, the earth is not flat, so the ground track is not a straight line. The situation is, in fact, better expressed by the diagram below.

For the purpose of this discussion, we will assume that the earth is a sphere of radius $R$. We will also assume that the satellite is in circular orbit around the earth, where the radius of that orbit is $r$. This means that:

$$ r = R + h $$ We will take this opportunity to define a ratio $b$ that will come in handy later: $$ b \equiv \frac {R} {r} $$

This diagram introduces the useful quantity central angle $\theta$ which is the angular distance between the observer and the satellite as seen from the center of the earth. For starters, as long as $\theta$ is expressed in radians, we can see that:

$$ D = R \cdot \theta $$ Also, some simple trigonometry leads us to the relation: $$ d = r \cdot \sqrt{1 + b^2 - 2b \cdot \cos \theta} $$

Alternatively, if we want to use the elevation $\beta$, taking a little help from the dotted extensions to the line segments in the figure above we can figure out that:

$$ d = r \cdot \left[ \sqrt{1 - b^2 \cdot \cos^2 \beta} - b \cdot \sin\beta \right] $$

Taking some more help from the dotted extensions to the line segments in the figure above we can deduce that:

$$ \cos {\Big( \beta + \theta \Big)} = b \cdot \cos \beta \qquad \implies \tan \beta = \frac {\cos \theta - b} {\sin \theta} $$ $$ r\cdot\sin\theta = d\cdot\cos\beta $$

It should be easy to see that as long as $\beta$ lies between $0^\circ$ and $+90^\circ$, $\theta$ will also be between $0^\circ$ and $+90^\circ$. In fact, we can put tighter bounds around $\theta$:

$$ 0 \le \theta \le \cos^{-1} b $$

But we still need find $\theta$ to actually compute $d$, $D$ and $\beta$.

It turns out that to obtain the values of $\theta$ and $\psi$ we need to use spherical trigonometry and great circle navigation from which we get the following results:

$$ \cos \theta = \sin \varphi \cdot \sin \varphi_s + \cos \varphi \cdot \cos \varphi_s \cdot \cos {\Delta\lambda} $$ $$ \tan \psi = \frac {\sin {\Delta \lambda} } {\cos \varphi \cdot \tan \varphi_s - \sin \varphi \cdot \cos {\Delta\lambda}} $$
Caveat: At this point we have to take note that for every valid value of a trigonometric ratio $\sin \phi$, $\cos \phi$ or $\tan \phi$, there are two valid values of the angle $\phi$ between $-180^\circ$ and $+180^\circ$, each in a different quadrant. If the ratio in question is positive, one of the values of $\phi$ lie in the first quadrant (between $0^\circ$ and $+90^\circ$). If we are sure that is the one we want we can safely use the $\mathrm{sin}^{-1}$, $\mathrm{cos}^{-1}$ or $\mathrm{tan}^{-1}$ function respectively to extract the value of $\phi$. In fact: $$ \sin^{-1} \sin \phi = \phi \implies -90^\circ \le \phi \le +90^\circ $$ $$ \cos^{-1} \cos \phi = \phi \implies 0^\circ \le \phi \le +180^\circ $$ $$ \tan^{-1} \tan \phi = \phi \implies -90^\circ \le \phi \le +90^\circ $$ However, when value of angle $\phi$ can legitimately be anything between $-180^\circ$ and $+180^\circ$, if we simplistically use the function $\mathrm{tan}^{-1}$ (which returns values between $-90^\circ$ and $+90^\circ$) we may end up with the wrong value. Therefore, in order to obtain the correct value of $\phi$, we have to back up a bit and use the two-argument version of $\mathrm{tan}^{-1}$, for which: $$ \tan^{-1} { \left[ \begin{matrix} {\color{red}{\sin \phi}} \\ {\color{blue}{\cos \phi}} \end{matrix} \right] } = \phi \implies -180^\circ \le \phi \le +180^\circ $$ This function is available as atan2 in several computer languages. The function $\tan^{-1} \left[ \begin{matrix} \color{red}{x} \\ \color{blue}{y} \end{matrix} \right]$ is $\texttt{atan2} \left( \color{red}{\texttt{x}} , \color{blue}{\texttt{y}} \right)$ in C, $\texttt{std::atan2} \left( \color{red}{\texttt{x}} , \color{blue}{\texttt{y}} \right)$ in C++, and $\texttt{Math.atan2} \left( \color{red}{\texttt{x}} , \color{blue}{\texttt{y}} \right)$ in Java. The arguments are reversed in Microsoft Excel and SQL, where it is $\texttt{ATAN2} \left( \color{blue}{\texttt{y}} , \color{red}{\texttt{x}} \right)$.

Since we have established that we are only interested in values of $\theta$ and $\beta$ that are between $0^\circ$ and $+90^\circ$, we can safely say:

$$ \theta = \cos^{-1} { \Big( \sin \varphi \cdot \sin \varphi_s + \cos \varphi \cdot \cos \varphi_s \cdot \cos {\Delta\lambda} \Big) } $$ $$ \beta = \tan^{-1} {\frac {\cos \theta - b} {\sin \theta}} $$ We can even deduce that $0^\circ \le \theta + \beta \le 90^\circ$, letting us safely say: $$ \theta = \cos^{-1} { \Big( b \cdot \cos \beta \Big) } - \beta $$

However, the value of angle $\psi$ can legitimately be anything between $-180^\circ$ and $+180^\circ$, so:

$$ \psi = \tan^{-1} { \left[ \begin{matrix} {\color{red}{\sin {\Delta \lambda} \cdot \cos \varphi_s}} \\ {\color{blue}{\cos \varphi \cdot \sin \varphi_s - \sin \varphi \cdot \cos \varphi_s \cdot \cos {\Delta\lambda}}} \end{matrix} \right] } $$

So far, we have referred to the object in the sky as a satellite, but even though we have mentioned its orbit we haven't actually used the orbit in any way. This makes all of the math above applicable not only to a satellite but also to a star, the sun, the moon, a meteor, a plane, a bird or Superman. That changes in the section where we consider ground tracks.

Looking down from the satellite

What about the satellite's point of view? The satellite will see the Earth as a disc of an angular diameter:

$$ 2\cdot \sin^{-1} b $$

The part of the Earth's surface that is visible from the satellite is the 'field of view' of the satellite, and it lies inside a circle that is drawn on the surface of the earth and is defined by all the points for which $\beta = 0^\circ$. The distance along the ground from the subsatellite point to each point on this circle is:

$$ R\cdot\cos^{-1} b $$

And the longest distance between the satellite and a point on earth visible from the satellite is to points on that circle:

$$ R \cdot \sqrt {\frac {1} {b^2} - 1} $$

The shortest distance between the satellite and a point on earth visible from the satellite is, of course, the distance $h$ to the subsatellite point, and that equals:

$$ R \cdot \left(\frac 1 b - 1 \right) $$

This 'field of view' circle covers the following fraction of the surface of the earth:

$$ \%_1 = \cfrac {1 - b} {2} $$

Orbits and the Ground Track

The path traced out on the earth by the subsatellite point as the satellite goes around its orbit around the Earth is the orbit's ground track.

We begin by defining the orbit to be a closed path around a spherical Earth, such that every point in that path and the center of the Earth lie on the same plane. The angle $i$ between this plane (the orbital plane) and the plane of the Earth's equator is the orbit's inclination. $0^\circ \le i \lt 90^\circ$ when the satellite orbits 'prograde' (in the same direction as the Earth's rotation), and $90^\circ \lt i \le 180^\circ$ when the satellite orbits 'retrograde' (in the opposite direction to the Earth's rotation.) When the orbit passes directly over the poles, we have $i = 90^\circ$.

First Approximation

As a first approximation, we assume that the Earth is not rotating. Or to be more pedantic, since we need the Earth's rotation to define its equator, we assume that it is rotating so slowly that its rotation can be ignored.

It should be easy to satisfy ourselves that the ground track will be the circle where the orbital plane and the spherical Earth intersect. Since the plane passes through the center of the Earth it divides the Earth into two equal hemispheres, and this in fact makes the ground track a great circle. What does the ground track of this orbit look like on a map? We can find that out by plotting the latitude $\varphi_s$ and longitude $\lambda_s$ of each possible subsatellite point.

In the diagram above, to the left we have a side view of the Earth. To the right, the figure at the top is the top view of the latitude $\varphi$, while the figure at the bottom is the corresponding top view of the satellite's ground track. Considering the Purple, Red & Green line segments, we see that the length of each one of them can be expressed in two different ways. Equating them, we get: $$ R \cdot \sin \gamma = R \cdot \cos \varphi \cdot \sin \lambda $$ $$ R \cdot \cos \gamma \cdot \cos i = R \cdot \cos \varphi \cdot \cos \lambda $$ $$ R \cdot \cos \gamma \cdot \sin i = R \cdot \sin \varphi $$ Combining the first two equations, we can immediately see that: $$ \tan \lambda = \frac {\tan \gamma} {\cos i} $$ Remembering the caveat about quadrants, we rephrase this slightly to obtain: $$ \lambda = \tan^{-1}{ \left[ \begin{matrix} \color{red}{\sin \gamma} \\ \color{blue}{\cos \gamma \cdot \cos i} \end{matrix} \right]} $$ $$ \varphi = \sin^{-1} {\Big( \cos \gamma \cdot \sin i \Big)} $$

We can plot $\varphi$ against $\lambda$ to get something like this picture (this one assumes $i=60^\circ$):

This would be the ground track if the Earth wasn't rotating. But since the Earth does in fact rotate, and the longitude reference meridian rotates with it, we also need to correct for Earth's rotation. Considering that the Earth rotates at a constant speed with a period $T_e$, we can compute a correction factor which is the amount $\delta$ the Earth rotates in the time $\tau$ it takes the satellite to cover the angle $\gamma$ in its orbit:

$$ \delta = \frac {2\pi} {T_e} \cdot \tau $$

For any given value of $\gamma$, then, the latitude $\varphi_s$ of the ground track point will be the $\varphi$ computed above while the longitude $\lambda_s$ will be $\lambda - \delta$.

Second Approximation: Circular Orbit around a Rotating Earth

If the orbit is a circular path that is concentric (i.e. shares the center) with the Earth, and the satellite orbits the Earth at constant speed with a period $T$, we have:

$$ \delta = \frac {2\pi} {T_e} \cdot \underbrace{\frac T {2\pi} \cdot \gamma}_{\tau} = \gamma \cdot \frac {T} {T_e} $$

Therefore, for one full orbit where $\gamma$ goes from $-180^\circ$ to $+180^\circ$, we will have:

$$ \varphi_s = \sin^{-1} {\Big( \cos \gamma \cdot \sin i \Big)} $$ $$ \lambda_s = \tan^{-1} { \left[\begin{matrix} \color{red}{\sin \gamma} \\ \color{blue}{\cos \gamma \cdot \cos i} \end{matrix}\right]} - \gamma \cdot \frac {T} {T_e} $$

What is the direction of the subsatellite point along ground track as the satellite goes around the orbit? As the satellite moves, $\gamma$ increases. A little calculus will tell us how $\lambda_s$ behaves as that happens:

$$ \frac {d \lambda_s} {d\gamma} = \frac {\sec i} {1 + \tan^2 i \cdot \sin^2 \gamma} - \frac {T} {T_e} = \frac {\cos i} {\cos^2 i + \sin^2 i \cdot \sin^2 \gamma} - \frac {T} {T_e} $$

Depending on whether this derivative is positive or negative the satellite appears to be moving East (prograde), or West (retrograde). For all $90^\circ \lt i \le 180^\circ$ this value is always negative, which does match our 'common sense' expectation that being East-to-West orbits to begin with, such satellites should always appear to move retrograde.

Also, when $i=90^\circ$ this derivative equals $- \cfrac T {T_e}$ which is always negative.

However, for all $0^\circ \le i \lt 90^\circ$, things are somewhat complicated. Since:

$$ \cos i \le \frac {\sec i} {1 + \tan^2 i \cdot \sin^2 \gamma} \le \sec i $$

We can deduce that:

$0^\circ \le i \lt 90^\circ$ $i = 90^\circ$ $90^\circ \lt i \le 180^\circ$
$\cfrac {T} {T_e} \lt \cos i$ $\cos i \le \cfrac {T} {T_e} \le \sec i$ $\cfrac {T} {T_e} \gt \sec i$
The satellite always appears to move prograde. The satellite appears to move prograde over high latitudes, but it appears to move retrograde close to the equator, specifically when
$\left| \varphi_s \right| \lt \sin^{-1} \sqrt {1 - \cfrac {T_e} {T} \cdot \cos i}$
The satellite always appears to move retrograde.

When we have a ground track that is sometimes prograde and some times retrograde, it is possible to have ground tracks that are 'loopy' like the figure below:

Final Approximation: Elliptical Orbit

In general an orbit is elliptical with the Earth at its focus. For the kind of stuff I'm discussing in this post, a satellite's orbit can be described by five quantities in general:

  1. The average distance $a$ between the satellite and the center of the Earth: this is the only orbital parameter on which the period $T$ of the orbit depends. Satellites orbits are in general elliptical, so this quantity is the semi-major axis of that ellipse. A circle is a special case of an ellipse where the two foci and center of the ellipse coincide: in that case this quantity maps to the radius of that circle.
  2. The inclination $0^\circ \le i \le 180^\circ$ of the orbit. When $i=0^\circ$ and $i=180^\circ$ we have equatorial orbits.
  3. The eccentricity $0 \le e \lt 1$ of the orbit. When $e=0$ we have circular orbits. The distance $r$ between the satellite and the center of the Earth is bounded $a\cdot\left(1 - e\right) \le r \le a\cdot\left(1 + e\right)$, which means that $a\cdot\left(1-e\right) \gt R$
  4. When the orbit is both elliptical ($0 \lt e \lt 1$) and inclined ($0^\circ \lt i \lt 180^\circ$) we also need to consider the the angle (as seen from the center of the Earth) between two points in that orbit:
    • One point that is a characteristic of the inclined nature of the orbit: we can pick any one of the northernmost point, the southernmost point, or any of the two points where the orbit crosses the equator (called the ascending node and descending node).
    • One point that is a characteristic of the elliptical nature of the orbit: we can choose either the point where the orbit comes closest to Earth (the perigee) or the point where it goes farthest away from Earth (the apogee).
    I choose the angle $A$ between the northernmost point of the orbit and the perigee because I've been using the northernmost point of the orbit as a reference point in this post. Conventionally, the angle $\omega$ between the ascending node and perigee is used for this, and this angle is 90° more than $A$.
  5. Finally, to completely specify the orbit we need to nail down the actual direction the orbit is inclined in. Being a direction in 3 dimensions, we need two angles, and the inclination is one of them. For the other, if I needed it for any derivation, in keeping with what we have been using so far I would use the right ascension of the northernmost point of the orbit $B$. Conventionally, the right ascension $\Omega$ of the ascending node is used for this, and it is 90° less than $B$.

Right ascension is a kind of celestial longitude: since the Earth rotates we cannot use the regular longitude to pinpoint things in space. The reference direction used is the line along which the Earth's equatorial plane and its orbital plane (around the sun) intersect. To be more precise, it is the direction that points from the Earth to the Sun at the Vernal Equinox.

With elliptical orbits $r$ varies, and the speed of the satellite along the orbit also varies. Since we have been counting $\gamma$ from $\left(i,0\right)$ the northernmost point in the orbit, we can first define $\alpha$:

$$ \alpha \equiv \gamma - A $$ Then, we have: $$ r = a\cdot \frac {1 - e^2}{1 + e \cdot \cos {\alpha} } $$ $$ \tau = \frac {T}{2\pi} \left( \tan^{-1} \left[\begin{matrix} \color{red}{ \sqrt {1 - e^2} \cdot \sin \alpha} \\ \color{blue}{e + \cos \alpha} \end{matrix}\right] - \frac {e \cdot \sqrt{1 - e^2} \cdot \sin \alpha} {1 + e \cdot \cos \alpha} \right) + \frac {T}{2\pi} \underbrace { \left( \tan^{-1} \left[\begin{matrix} \color{red}{ \sqrt {1 - e^2} \cdot \sin A} \\ \color{blue}{e + \cos A} \end{matrix}\right] - \frac {e \cdot \sqrt{1 - e^2} \cdot \sin A} {1 + e \cdot \cos A} \right) }_{\equiv \; \mathcal{C},\; \textsf{a constant for this orbit}} $$

Clearly, the rate of change of $\gamma$ with time $\tau$ is no longer the constant $\frac {2\pi} T$. Putting it all together, we get the ground track (for one full orbit where $\gamma$ goes from $-180^\circ$ to $+180^\circ$) as follows:

$$ \varphi_s = \sin^{-1} {\Big( \cos \gamma \cdot \sin i \Big)} $$ $$ \lambda_s = \tan^{-1} { \left[\begin{matrix} \color{red}{\sin \gamma} \\ \color{blue}{\cos \gamma \cdot \cos i} \end{matrix}\right]} - \frac {T}{T_e} \cdot \tan^{-1} \left[\begin{matrix} \color{red}{\sqrt {1 - e^2} \cdot \sin {\left(\gamma - A\right)} } \\ \color{blue}{e + \cos {\left(\gamma - A\right)}} \end{matrix}\right] + \frac {T}{T_e} \cdot \frac {e \cdot \sqrt{1 - e^2} \cdot \sin {\left(\gamma - A\right)}} {1 + e \cdot \cos {\left(\gamma - A\right)}} - \frac {T} {T_e}\cdot\mathcal{C} $$

This is as far as we can go with trigonometry and the understanding of orbits up to Kepler (1571-1630). The only bit of physics required to explain what we've discussed is a 'central force' that is inversely proportional to the square of the distance from the focus and is pulling the satellite towards the focus.

Predictable Complications

The explanation of what the 'central force' really was came from Newton (1642-1727). It turned out to be, as we all know by now, gravitation. However, since the Earth is neither a point mass nor a perfect sphere of uniform density, the force of gravitation experienced by the satellite does not always point towards the center of the Earth. The direction it does point to depends on the shape of the Earth, or more precisely the shape its gravitational field. Even more precisely the direction of the force at any point is perpendicular to the equipotential surface (the potential in question being the gravitational potential) passing through that point. So it is the shape of this equipotential surface we are after. What is this shape?

  • The first approximation is a sphere.
  • The second approximation is an 'oblate spheroid' (a.k.a. an 'ellipsoid of revolution') where the radius at the equator is a little bit more than that at the poles.
  • The third approximation is a 'tri-axial ellipsoid', where even the equator is somewhat elliptical, and thus the radii in the $x$, $y$ and $z$ directions are all somewhat diferent.

... but of course these are all approximations. The shape is officially a 'Geoid' which simply is another way of saying 'earth-shaped' and thus doesn't tell us anything we don't know. Crucially, it isn't the shape of the actual surface of the Earth. Instead, imagine if the Earth could be entirely covered, up to the tip of Mt. Everest, by a liquid of negligible density: the shape of the Geoid would then be be the shape of the 'sea level' of this 'ocean'.

The shape of the Geoid is, of course, just what it happens to be: we can only determine it by taking a lot of measurements. This has been done, and continues to be done. But of course it is just a large volume of data, and we can try to make it tractable by trying to reduce to a set of parameters of a spherical harmonic series. Here is the series for the Earth's gravitational potential $U$ at any point $\left(r, \varphi, \lambda \right)$, where each $P_n $ is a Legendre Polynomial and each $P^m_n$ is an Associated Legendre Polynomial:

$$ U \left( r, \varphi, \lambda \right) = \frac {G\cdot M_e} r + \sum_{n=2}^{\infty} { \cfrac { J_{n} \cdot P_n \left( \sin \varphi \right) } { r^{n+1} } } + \sum_{n=2}^{\infty} { \sum_{m=1}^{n} { \cfrac {J_{nm} \cdot P^m_n \left( \sin \varphi \right) \cdot \cos {\left(m\cdot\left[ \lambda - \lambda_{nm}\right]\right)}} { r^{n+1} } } } $$

For a perfect sphere, the $J_n$ and $J_{nm}$ terms would all be zero. But for the real Earth, they aren't. The coefficient with the largest impact is $J_2$, and it quantifies the 'oblate spheroid' shape. The effects of $J_2$ are:

  • The slow rotation of the orbital plane itself around the Earth's axis (Nodal Precession). For every turn of the orbit, this causes a shift in the ascending node (and descending node) by: $$ \Delta\Omega_{J_2} = -3\pi\cdot \cfrac {J_2} {G\cdot M_e} \cdot \left[\frac {1} {a \cdot \left(1-e^2\right)}\right]^2 \cdot \cos i $$
  • The slow rotation of the major axis (line joining the points on the orbit closest and farthest from the Earth) in the orbital plane (Apsidal Precession). For every turn of the orbit, this causes a shift in the apogee (and perigee) by: $$ \Delta\omega_{J_2} = -6\pi\cdot \cfrac {J_2} {G\cdot M_e} \cdot \left[\frac {1} {a \cdot \left(1-e^2\right)}\right]^2 \cdot \left(\frac 5 4 \cdot \sin^2 i - 1 \right) $$

The value of $\Delta\Omega_{J_2}$ becomes zero for $i = 90^\circ$ (polar orbits), while $\Delta\omega_{J_2}$ becomes zero for $i = \tan^{-1} 2 \approx 63.435^\circ$. Other than that, however, both quantities increase as the orbit becomes more eccentric, and as the average height of the orbit reduces.

The next coefficient to have an impact is $J_{22}$ (and its associated longitude $\lambda_{22}$). It quantifies the very slight ellipticity of the Earth's equator. As it turns out, for the Earth $\lambda_{22}\approx-15^\circ$, so the qualitative impact of the parameter $J_{22}$ is a very gentle nudge towards longitudes $+75^\circ/-105^\circ$, which are points of stability $90^\circ$ East and West of $-15^\circ$.

Finally, Einstein (1879-1955) also makes an appearance here. There is also some Apsidal Precession entirely independent of the shape of the Earth, and that is due to General Relativity. This is the following for every turn of the orbit: $$ \Delta\omega_{GR} = +2\pi \left[ \frac {\Huge{1}} {\sqrt{1 - 6\cdot \cfrac {G\cdot M_e} {c^2} \cdot \cfrac 1 {a\cdot\left(1 - e^2\right)} } } - 1\right] $$ ... which is approximately: $$ +6\pi\cdot \frac {G \cdot M_e} { c^2 } \cdot \left[\frac {1} {a \cdot \left(1-e^2\right)}\right] $$ ... when $a\left(1-e^2\right) \gg \frac {G\cdot M_e}{c^2}$ which it usually is, making the effect negligible for all practical purposes.

Orbits with Fixed Ground Tracks

Now let us consider orbits where the satellite passes only over a fixed track on the surface of the Earth

Obviously, any orbit where the orbital plane is the same as the Earth's equatorial plane ($i=0^\circ$) satisfies this criterion, and their ground track is the Equator itself: this is an equatorial orbit. But this is also possible with an inclined orbit, provided the time period of the orbit is a rational multiple of the length of the Earth's period of rotation. That is to say, when $p$ and $q$ are positive integers with no common factors, the time period $T$ of the orbit is:

$$ T = {p \over q} \cdot T_e $$

When $p \neq q$, we can have the following situations:

$0^\circ \le i \lt 90^\circ$
(Prograde)
$q \gt p$
(Subsynchronous)
The orbit will repeat its ground track on the map in groups of $q - p$ passes, each from longitude $-180^\circ$ to longitude $+180^\circ$. For every $p$ rotations of the Earth the satellite goes around it $q$ times, but as seen from the Earth it appears to have gone around $q - p$ times in the prograde (W → E) direction.
$p \gt q$
(Supersynchronous)
The orbit will repeat its ground track on the map in groups of $p - q$ passes, each from longitude $+180^\circ$ to longitude $-180^\circ$. For every $p$ rotations of the Earth the satellite goes around it $q$ times, but as seen from the Earth it appears to have gone around $p - q$ times in the retrograde (E → W) direction.
$90^\circ \le i \lt 180^\circ$
(Retrograde)
The orbit will repeat its ground track on the map in groups of $q + p$ passes, each from longitude $+180^\circ$ to longitude $-180^\circ$. For every $p$ rotations of the Earth the satellite goes around it $q$ times, but as seen from the Earth it appears to have gone around $q + p$ times in the retrograde (E → W) direction.

In the picture below, we show some such prograde orbits, all inclined $i=60^\circ$ to the equatorial plane, where $\left|p - q\right| = 1$.

In fact, the satellites of the GPS constellation are in orbits where $T = \frac 1 2 T_e$, except that their inclination is a little less ($55^\circ$).

When $p = q$ for a prograde orbit we have $T = T_e$: such orbits are called synchronous orbits. As we can see in the following figure, the ground track of a synchronous orbit is bounded in extent within the Earth's latitude/longitude grid, staying between latitudes $+i$ and $-i$ and between longitudes that are determined by the orbit's parameters. When inclined these orbits have a 'figure of eight' pattern fixed over the Earth. When some eccentricity is sometimes added to such orbits it slightly warps the 'figure of eight' shape towards some direction.

Inclined synchronous orbits are a favorite with regional satellite positioning systems, such as the Indian NAVIC/IRNSS ($i=29^\circ$) and the Chinese Beidou ($i=55^\circ$) systems.

Synchronous orbits that are elliptical but not inclined (i.e. $0 \lt e \lt 1$ and $i=0^\circ$) have ground tracks that are arcs of the equator: the satellite moves to and fro, East to West and back, over such arcs.

When we have a synchronous orbit which is circular ($e=0$) and equatorial ($i=0^\circ$), the ground track is a single point on the equator: the satellite appears to 'stand still in the sky' over this point, and consequently this orbit is the stationary orbit.

Circular Equatorial Orbits

An equatorial orbit is one with an inclination of $0^\circ$. Since all circular orbits around the earth share the center with the earth itself, this puts all points of such orbit directly above the equator, and sure enough substituting $i=0$ into $\sin^{-1} \left( \cos \gamma \cdot \sin i \right)$ does confirm that the latitude of each point on the ground track of such an orbit is $0^\circ$.

Since $\varphi_s = 0$, the key numbers for such orbits are as follows:

$$ \psi = \tan^{-1} {\left[ \begin{matrix} \color{red}{\sin{\Delta\lambda}} \\ \color{blue}{ -\sin \varphi \cdot \cos {\Delta\lambda}} \end{matrix}\right]} $$ $$ \beta = \tan^{-1} \frac {\cos \varphi \cdot \cos {\Delta\lambda} - b} {\sqrt{1 - \cos^2 \varphi \cdot \cos^2 {\Delta\lambda}}} $$ $$ d = r \cdot \sqrt{1 + b^2 - 2b \cdot \cos \varphi \cdot \cos {\Delta\lambda} } $$ $$ D = R \cdot \cos^{-1} {\Big( \cos \varphi \cdot \cos {\Delta\lambda} \Big)} $$

For any satellite in such an orbit to be visible by an observer, the following conditions must be satisfied. First, the observer must be located in a low enough latitude:

$$ \left| \varphi \right| \le \cos^{-1} b $$ Next, the satellite's longitude must not be too far East or West of the observer's: $$ \left| \Delta\lambda \right| \le \cos^{-1} {\frac {b} {\cos \varphi} } $$ When the conditions above are met, the satellite can be seen no higher in the sky than: $$ 0 \le \beta \le \tan^{-1} {\frac {\cos \varphi - b} {\sin {\left| \varphi \right|} } } $$ Finally, we consider the slice of horizon $\Delta\psi$ outside of which the satellite cannot be seen. $$ \Delta\psi = 2 \cdot \tan^{-1} \frac {\sqrt{\cos^2 \varphi - b^2}} {b \cdot \sin {\left| \varphi \right|} } $$

The meaning of $\Delta\psi$ is explained in the figure below. This is the difference in azimuth between the two points on the horizon for which $\beta = 0$, and $0 \le \Delta\psi \le 180^\circ$.

From the point of view of a satellite in a circular equatorial orbit, the following fraction of the Earth's surface will always be invisible:

$$ \%_2 = 1 - \sqrt{1 - b^2} $$

A matter of timing

For any object (such as a satellite) the mass $m$ of which is much smaller than the mass of the earth $M_e$, the orbital period $T$ is given by:

$$ T = 2 \pi \cdot \sqrt { \frac {a^3} {G \cdot M_e} } $$

As long as we are interested only in circular orbits (where $r = a$) and if we consider the Earth to be a sphere with uniform density $\rho$, we can express this in terms of $b$:

$$ M_e = \rho \cdot \frac {4\pi} {3} \cdot R^3 \implies T = \sqrt { \frac {3\pi} {G \cdot \rho \cdot b^3} } $$

It is tempting to try to define a 'visibility time' $t$ which is defined, for an observer located on the ground track of the satellite, as the time between the 'rising' and 'setting' of the satellite. If we assume a non-rotating Earth, this figure is quite easy to calculate, being the twice time taken by the satellite to cover twice the central angle between the $\beta = 0^\circ$ (horizon) and $\beta = 90^\circ$ (zenith) positions from the point of view of the observer:

$$ t = 2 \cdot \cos^{-1} b \cdot \frac T {2\pi}= \sqrt{\frac 3 {\pi \cdot G \cdot \rho} } \cdot \sqrt{ \frac 1 {b^3} } \cdot \cos^{-1} b $$

The Earth does rotate, and once we consider the complicated ground tracks this leads to for inclined orbits we realize that in general the value of $t$ will be different for different points on the ground track. However, if we can confine ourselves to equatorial circular orbits we can obtain values of $t$ that are the same at every point on the ground track:

Prograde equatorial orbit
($i = 0^\circ$)
Retrograde equatorial orbit
($i = 180^\circ$)
$$ t = \frac {\cos^{-1} b} { \left| \sqrt {\cfrac {\pi \cdot G \cdot \rho} {3}} \cdot \sqrt {b^3} - \cfrac {\pi} {T_e} \right| } $$ $$ t = \frac {\cos^{-1} b} { \sqrt {\cfrac {\pi \cdot G \cdot \rho} {3}} \cdot \sqrt {b^3} + \cfrac {\pi} {T_e} } $$

Putting this in a graph (spoiler alert: we've used ${T_e}^2\cdot G \cdot \rho = 2732$ which is actually the value for our Earth) we get something like this where the red line represents a retrograde orbit while the blue one stands for a prograde orbit:

That part of the blue line where $t$ goes off to infinity corresponds to when $b$ is:

$$ \sqrt[\leftroot{3}\uproot{4}\scriptstyle 3] {\cfrac {3\pi} {{T_e}^2\cdot G \cdot \rho}} $$

... and that is going to the topic of the following section

The Geostationary Orbit

The geostationary orbit is the Earth's stationary orbit. It is a prograde (West-to-East) circular equatorial orbit the time period of which is exactly equal to the Earth's period of rotation, also known as the length of a sidereal day: 23 hours, 56 minutes, 4.09 seconds. In other words, it shares the center, axis and angular velocity with the earth itself. As a consequence, from the surface of earth any satellite in this orbit appears to be stationary at one point in the sky: it doesn't rise, it doesn't set, and so $t = \infty$.

The relevant parameters of our Earth are as follows: $$ R = 6371 \,\mathrm{km}, \quad T_e = 23.93447 \,\mathrm{hr} = 86164.09\,\mathrm{s}, \\ \quad G \cdot M_e = 3.986 \times 10^{14} \;\mathrm{m^3 s^{-2}}, \quad G \cdot \rho = 3.68 \times 10^{-7} \;\mathrm{s^{-2}} $$
The value of $G \cdot M_e$ and $G\cdot\rho$ can be deduced from the more commonly known value of $9.81\,\mathrm{m\,s^{-2}}$ for $g$, the standard acceleration due to gravity, because: $$ g = \frac {G \cdot M_e} {R^2} $$

From this information, we can deduce the value of $b$ for the geostationary orbit, and from that its radius too:

$$ b = 0.1511, \quad r = 42164 \;\mathrm{km} $$
This implies that $h = 42164 - 6371 = 35793 \,\mathrm{km}$. However, the geostationary orbit is commonly referred to as one that is $35786\,\mathrm{km}$ above the earth. This discrepancy is because while in this article we have considered the earth to be spherical, it is actually better described as an ellipsoid, with a polar radius of $6357\,\mathrm{km}$ and an equatorial radius of $6378\,\mathrm{km}$. Since the geostationary orbit is entirely over the equator, it is indeed $42164 - 6378 = 35786 \,\mathrm{km}$ over the earth.

Since $\cos^{-1} {0.1511} = 81.3093^\circ$, we can infer that the geostationary orbit will only be visible to observers between latitudes 81.3093°N and 81.3093°S.

The 'field of view' as seen from the satellite is a circle, every point in the circumference of which is $9041\,\mathrm{km}$ away along the ground from the subsatellite point: this circle contains all points for which $\beta = 0^\circ$, and covers 42.4% of the surface of the Earth. However, 1.1% of the Earth, half of it around the North Pole and the other half around the South Pole, will always be outside the field of view of any geostationary satellite.

If we draw such circles for progressively higher values of $\beta$, they will form a bulls-eye pattern centered on the subsatellite point. When we draw that pattern on a flat map using a cylindrical projection, we get something like this:

This means, of course, that inside of those circles, the look angle to the satellite is at least the corresponding value for $\beta$.

$\beta$ $\left| \varphi \right|$ $D$ $d$ $\%_{1}$ $\%_{2}$
81.3° 9041 km 41680 km 42.4% 1.1%
15° 66.6° 7406 km 40064 km 30.1% 8.2%
30° 52.5° 5836 km 38616 km 19.5% 20.7%
45° 38.9° 4322 km 37418 km 11.1% 37.2%
60° 25.7° 2854 km 36526 km 4.9% 56.7%
75° 12.8° 1419 km 35978 km 1.2% 77.9%

The numbers corresponding to the bulls-eye pattern, going from the outside inwards to the center, are as in the table above.

Satellites and the Consumer

Antennas, Frequency Bands and Orbit Types

Small Fixed Dish, Omnidirectional, Flat-panel Phased Array (Large Dish, Gimbal-mounted Dish)

VHF, L, S, C, Ku, Ka, V

LEO/MEO, GEO, IGSO (HEO both meanings)

Direct Broadcast

Direct broadcast is a form of one-way communication. The typical use is to broadcast 'linear' audio and video programming (old fashioned TV & radio, only over satellite.) It is a very good fit for GEO:

  1. Latency is not a problem
  2. The end user does not transmit, and is thus not bothered by the free space path loss over the 36000+ km to the satellite
  3. Relatively few strategically placed GEO satellites can easily cover pretty much every place on earth where people live.

When we consider GEO satellites, for broadcast receivers in homes the spatial relationship between the receive antenna and the satellite never changes, so small fixed dish antennas are the best choice. For receivers in vehicles the azimuth (and to a lesser extent the elevation) to the satellite changes as the vehicle moves, so the upcoming flat panel phased arrays will be necessary.

Navigation/Positioning

Satellite navigation is in fact a form of broadcast: each navigation satellite repeatedly broadcasts the following:

  1. The current time $t_1$ read off a highly precise atomic clock on board the satellite.
  2. The orbital parameters of the satellite itself (the equivalent of the four we discussed earlier, and a few more) at a reference time $t_0$ not too long ago in the past.

A receiver, of course, receives this message at a later time $t_1+\Delta t$, where $\Delta t$ is time taken by radio waves (which travel 300 kilometers every millisecond) to travel from the satellite to the receiver and is unknown at the moment of reception. The receiver can, however, use this data to compute the position of the satellite at the moment $t_1$ when this message was transmitted. After it has such gathered such data from four different satellites, it can solve for its own exact position in 3 space dimensions $\left(\varphi, \lambda, h\right)$ as well as the precise current time $t$.

System MEO IGSO GEO Band
$i$ $h$
GLONASS Russia 64°8′ 19140 km - - L
GPS USA 55° 20180 km - - L
Beidou/COMPASS China 55° 21500 km 55° Yes L
Galileo EU 56° 23222 km - - L
IRNSS/NAVIC India - - 29° Yes L, S

Compared to audio and video broadcast, however, the coverage requirements are much more stringent and the data rates/volumes are much lower. The application makes it impractical to use receiver antennas that require 'pointing' towards the satellite, so receivers use 'built-in' or external omnidirectional antennas. L-band and S-band frequencies make this possible while keeping the antenna size reasonable.

Two-way Communications

Most of today's bidirectional satellite communications take place over GEO satellites: a majority of the 400-plus satellites in the geostationary orbit today are used for some combination of bidirectional data communication and audio/video broadcast. However, there are a few challenges that make this less than ideal:

  1. The relatively long distance between the terminal and the satellite introduces a delay of between 120 & 140 milliseconds one way. In addition to the disconcerting effect this has on voice calls, this also degrades high bit-rate data streams that use acknowledgements.
  2. The distance also means a large value of Free Space Path Loss, which increases with the square of distance. Everything else being equal, this increases the uplink transmit power required at the terminal.
  3. In higher latitudes the look angle does get small, which means that in mountainous landscapes and downtown skylines the satellite may not be visible.

Some systems in non-geostationary orbits have also been deployed to address these shortcomings. However, since these satellites do not 'stand still in the sky' it is necessary to have a constellation of satellites to maintain continuous coverage as each individual satellite moves in and out of sight.

System $h$ $i$ Band
LEO Iridium 781 km 86.4° L
Orbcomm 825 km 45.0° VHF
Globalstar 1414 km 52.0° L
Gonets ~1500 km 82.5° VHF
MEO O3B 8063 km 0° (Equatorial) Ka

Other than the ones above, several others have been proposed:

System Status LEO/MEO IGSO GEO Band
$i$ $h$ $e$ $i$ $e$
ICO Abandoned 45.0° 10355 km - - - - L
Spaceway Incomplete - - - - - Yes Ka
55.0° 10352 km - - - -
Teledesic Abandoned 84.7° 1380 km - - - - Ka
Skybridge Abandoned 53.0° 1469 km - - - - Ku
Viasat Proposed 87.0° 8200 km - - - - Ka, V
O3B Proposed 0° (Equatorial) 8063 km - - - - Ka, V
70.0° 8063 km - - - -
Leosat Proposed 90° (Polar) 1400 km - - - - Ka
Kepler Proposed 90° (Polar) 600 km - - - - Ku
Oneweb (1) Proposed 87.9° 1200 km - - - - Ku
Oneweb (2) Proposed 87.9° 1200 km - - - - V
45° ± 1° 8500 ± 100 km - - - -
Boeing (1) Proposed 54.0° 1056 km - - - - V
- - - 63.4° (Critical) 0.2 -
Boeing (2) Proposed 45.0° 1030 km - - - - V
55.0° 1082 km - - - -
Telesat Proposed 99.5° 1000 km - - - - Ka, V
37.4° 1248 km - - - -
Spacex Proposed 53.8° 1110 km - - - - Ku
74.0° 1130 km - - - -
81.0° 1275 km - - - -
70.0° 1325 km - - - -
New Spectrum Proposed 63.4° (Critical) 13920 km 0.6 - - - C, Ku, Ka