# Tiangong-1 reentry impact place probabilities

Following all the news concerning Tiangong-1, I have seen the following plot made by the Aerospace Corporation where

• Blue area: null probability of impact
• Yellow area: high probability of impact
• Green area: low probability of impact

I wonder: why the yellow areas, that are more or less coincident when the space station flies over the maximum and minimum latitudes (see the groundtrack), are the most likely places to reentry?

• A first hint is the density of the lines depicting the flight path. Compare densities for Spain with the Galapagos Islands (near the -90 mark on the Equator). Three paths cross Spain, only one the Galapagos Islands. Second hint: note that due to the map projection, the yellow part is even stretched in the east/west direction. Actual path density is thus even larger in the yellow region. – Jens Mar 24 '18 at 11:49

From the pros, from the Spaceflight 101 article Tiangong-1 Re-Entry, click for full size: For a circular LEO, for $x_p$, $y_p$ in the plane of the orbit we can just write

$$x_p = \cos(\omega t)$$ $$y_p = \sin(\omega t),$$

and if it is inclined to the equator by an angle $i$, the $x$, $y$, $z$ coordinates when the $\hat{z}$ axis is parallel to the Earth's rotational axis will be

$$x = \cos(\omega t)$$ $$y = \sin(\omega t) \cos(i),$$ $$z = \sin(\omega t) \sin(i),$$

and the latitude will be

$$\lambda=\arctan\left(\frac{z}{\sqrt{x^2+y^2}} \right)$$

Here's a plot for inclinations of 20, 40, 60, and 80 degrees. Two Three things stand out.

1. As inclination approaches 90 degrees and the orbit becomes polar, the plot of latitude versus time becomes more triangular. Of course at exactly 90 degrees the latitude increases or degreases purely linearly with time for a circular orbit.
2. More interesting is that the U-shape of the time-binned histogram of latitude flattens out as well. For high inclination orbits, the amount of time spent per degree of latitude becomes much more even, except for the "ears" at max and min, where there it still sort-of "stalls" at the extrema before turning around again.
3. Even more interesting is the time-binned histogram rescaled by $1/ \cos(\lambda)$ for surface area rather than latitude, as recommended by @Litho's comment. If you were looking for debris, or looking to avoid getting hit by debris personally, this would be the plot for you. import numpy as np
import matplotlib.pyplot as plt

halfpi, pi, twopi = [f*np.pi for f in 0.5, 1, 2]

omega = twopi
N     = 20000
t     = np.linspace(0, 1, 20000)
incs  = [rads*d for d in (20, 40, 60, 80)]

lats = []
for inc in incs:
x = np.cos(omega*t)
y = np.sin(omega*t) * np.cos(inc)
z = np.sin(omega*t) * np.sin(inc)

rxy = np.sqrt(x**2 + y**2)
lat = np.arctan2(z, rxy)
lats.append(lat)

bins = np.arange(-90, 91, 1)

hists = []
for lat in lats:
latdegs = degs*lat
hists.append(np.histogram(latdegs, bins))

if True:
fig = plt.figure()

for lat in lats:
ax1.plot(t, degs*lat)
ax1.set_ylim(-90, 90)
ax1.set_title('latitude vs time', fontsize=16)

for a, b in hists:
ax2.plot(b[1:], a)
ax2.tick_params(labelleft='off')
ax2.set_xlim(-90, 90)
ax2.set_title('per unit latitude, area=1', fontsize=16)

for a, b in hists:
ax3.tick_params(labelleft='off')
ax3.set_xlim(-90, 90)
ax3.set_title('per unit area for @Litho, area=1', fontsize=16)

plt.show()

• Your histogram shows the time which a spacecraft spends over a given latitude. But what you need is the time it spends over a given unit area of surface over the course of many orbits, and to get it, you need to divide your histogram by the cosine of the latitude, since the area of Earth's surface between latitudes $\theta$ and $\theta + d\theta$, where $d\theta$ is very small, is approximately $2\pi R^2\cos(\theta)d\theta$, where $R$ is Earth's radius. – Litho Mar 23 '18 at 17:26
• I do not understand well the coordinate change you made. It should not be $x=\cos(\omega t)\cos(i)$ and $z=\sin(i)$? With $y$ coordinate I have no issues. – Julio Mar 23 '18 at 22:47
• @Julio I've plotted a circular orbit with the ascending node at $t=0$ and $x, y, z = 1, 0, 0$. At $t=0.25$, those will be $x, y, z = 0, \cos(i), \sin(i)$. All three must be time-dependent; $z=\sin(i)$ only would be unphysical. – uhoh Mar 24 '18 at 1:19
• @Litho I've explained the principle of the higher probability at the extrema, which is what the OP has actually asked about. Indeed if you look closely, those yellow bands are squiggly and represent an orbit propagation, but I'm not going to propagate a TLE of a tumbling, re-entering spacecraft weeks into the future and call it "correct". If altitude loss rate is uncertain (which it is), rotating-Earth-projected longitudes will also be smeared out with uncertainty since mean motion changes. My plot is correct as advertised but I've added a third plot for "per unit area" for you. See item #3. – uhoh Mar 24 '18 at 1:44

Simply put, that's where it spends most of the time. I think the chart is a bit exaggerated, but a satellite will be at it's northern and southern extrema far more. The pattern follows roughly a sinusoidal pattern, and is further enhanced because the circles shrink in size the further you go from the equator. I can't find a plot with the orbit, but here is a CDF of just the sinusoidal component. Keep in mind that this is further increased by the spherical nature of the globe. 