# Flying in Formation: The Orbital Dynamics of LISA’s Three Spacecraft

### A.  Introduction

The era of gravitational wave (GW) observational astronomy began on September 14, 2015, when LIGO’s (Laser Interferometer Gravitational Wave Observatory) two giant Michelson-like interferometers detected the gravitational radiation emitted as two black holes, with total mass of $\boldsymbol{70\, M_{\textup{Sun}}}$, collided and coalesced more than 1 billion light years from Earth. [Ref. 1]  This landmark discovery was followed by six others by LIGO and its sister interferometer VIRGO, located near Pisa, Italy.  In each case, the Earth-bound detectors captured a short burst of gravitational radiation emitted in the final seconds before two astronomical bodies — either black holes or neutron stars — collided and merged.  In those fleeting moments, the GW power radiated by the two colliding bodies was comparable to the light emitted by the entire visible universe!

Although GW astronomy is still in its infancy, scientists are already planning the next generation of GW observatories.  The huge LIGO detectors, with interferometer arms stretching for 4 km, are exquisitely sensitive, yet they are incapable of detecting other important astronomical events such as the merger of massive black holes, those with mass  $\boldsymbol{> 10^{4}\, M_{\mathrm{Sun}}}$.  The highest frequency of a GW emitted during the merger of two bodies is inversely proportional to the bodies’ total mass, and LIGO’s low frequency sensitivity is limited by its size and residual coupling to the low-pitched rumble of Earth’s seismic activity.  This is where LISA comes in.  Anticipated to launch in 2034, LISA (Laser Interferometer Space Antenna) is a much larger space-based GW observatory employing three identical spacecraft equipped with lasers, mirrors, telescopes and detectors to form a giant interferometer.  The spacecraft will be deployed about 0.35 AU from Earth and will fly in an equilateral triangular formation of side length $\boldsymbol{2.5\times 10^{6}\: \textup{km}}$ (Figure 1. )  Due to its enormous size and immunity to Earth’s vibrations, LISA will be sensitive to GW radiation at much lower frequencies than LIGO or VIRGO, thus complementing the Earth-bound instruments. [Ref. 2]

LISA is expected to gather data for 2-10 years.  To function for this length of time, its three spacecraft must maintain their separation and configuration with minimal assistance from on-board ion thrusters.  But simply injecting them into orbit with the desired spacing and orientation does not guarantee that they will remain so.  After injection, each spacecraft will move independently in its own elliptical orbit about the Sun, and unless the three orbits are chosen carefully, the configuration of the spacecraft will change unacceptably with time.

As an extreme example, consider the fleet of four spacecraft shown in Figure 2.  The “horizontal” and “vertical” spacecraft pairs represent two orthogonal arms of a GW interferometer, which are initially of length $\boldsymbol{L}$.  The four bodies lie in a plane that is inclined by angle φ to the plane of the ecliptic, and the center of the fleet moves about the Sun in a circular orbit of radius $\boldsymbol{R=1\: \textup{AU}}$.  Let $\boldsymbol{\phi =90^{\circ}}$ and assume $\boldsymbol{L\ll R}$.  (For LISA, $\boldsymbol{L/R=0.017}$.)  The distance of each spacecraft to the Sun is $\boldsymbol{r=\sqrt{R^{2}+L^{2}/4}\simeq R(1+L^{2}/8R^{2})}$, so to first order in $\boldsymbol{L/R}$$\boldsymbol{r=R}$ for each spacecraft.  Naïvely, it would seem that if all four spacecraft were to travel in circular orbits with period $\boldsymbol{T=1\: \textup{yr}}$, their configuration would be stable.  Clearly, if spacecraft C and D travel within the ecliptic plane, they would remain separated by $\boldsymbol{L}$ as they circle the Sun.  But the orbits of A and B are tilted from the ecliptic by the small angle $\boldsymbol{\delta =\pm L/2R}$, and after time $\boldsymbol{T/4=3\: \textup{mo}}$, they intersect and — in the worst case — the two spacecraft would collide!  A more physics-informed strategy is needed for a successful space-based interferometer. [Ref. 3]

In the following sections, we examine how the spacing and orientation of a fleet of spacecraft evolve with time, to understand the strategy underlying LISA’s design.  We begin by presenting an approximate treatment of motion in elliptical orbits, to derive expressions for the position of a spacecraft and the separation between spacecraft, as functions of time.  Using these expressions, we then carry out a straightforward numerical calculation of spacecraft position and separation vs. time, which will elucidate the orbital design adopted by LISA.  All of this can be accomplished using the elementary mathematics and physical principles found in the typical introductory calculus-level physics course.  Because of this and the current high interest in gravitational radiation, LISA offers instructors an ideal opportunity to illustrate basic orbital mechanics in an intriguing way.  The purpose of this post is to show one way to present this topic to undergraduates.

### B.  Position vs. Time in Elliptical Orbits

Each of LISA’s three spacecraft will revolve about the Sun in an elliptical orbit with period $\boldsymbol{T=1}\: \textup{yr}$, so each orbit must have a semi-major axis $\boldsymbol{a=1\: \textup{AU}}$.  If an orbit’s eccentricity $\boldsymbol{\epsilon \ll 1}$, as is the case for LISA, it can be approximated as a circle with the Sun offset from its center by a distance $\boldsymbol{\epsilon R}$, where $\boldsymbol{R}$ is the radius of the circle and also the semi-major axis of the ellipse ($\boldsymbol{R=a}$) (Figure 3).  [Ref. 4]

The position $\boldsymbol{\vec{r}(t)}$ of the orbiting body, relative to the Sun, is specified by its distance to the Sun $\boldsymbol{r}$ and the angle ν, called the true anomaly in astronomical literature.  Similarly, the position $\boldsymbol{\vec{R}(t)}$ of the body relative to the center of the circle is given by $\boldsymbol{R}$ and the angle θ, called the eccentric anomaly in the literature.  Let $\boldsymbol{t=0}$ and $\boldsymbol{\nu =\theta =0}$ when the body is at perihelion, its closest approach to the Sun.  The area swept out by $\boldsymbol{\vec{r}(t)}$ from perihelion to time $\boldsymbol{t}$ is equal to the circular sector swept out by $\boldsymbol{\vec{R}(t)}$ minus the shaded triangular area shown in Fig. 3.  Since the total area swept out by $\boldsymbol{\vec{r}(t)}$ in a full orbital period is $\boldsymbol{\pi R^{2}}$, Kepler’s 2nd (equal area) Law yields

$\LARGE\boldsymbol{\frac{t}{T}=\frac{\frac{1}{2}R^{2}\theta -\frac{1}{2}\epsilon R^{2}\mathrm{sin}\theta }{\pi R^{^{2}}}}$

or

$\LARGE \boldsymbol{\frac{2\pi t}{T}=\theta -\epsilon\, \mathrm{sin}\theta\equiv M}$,

where $\boldsymbol{M}$ is often called the mean anomaly.  (Note that if $\boldsymbol{\epsilon =0}$$\boldsymbol{M=\theta }$, so $\boldsymbol{M}$ is the angular position of a body moving in a circular orbit with period $\boldsymbol{T}$.)  This is a transcendential equation for $\boldsymbol{\theta (t)}$, and an approximate solution is easily found for $\boldsymbol{\epsilon \ll 1}$.  To first order in $\boldsymbol{\epsilon}$$\boldsymbol{\theta =M}$, so

$\large \boldsymbol{\theta(t) = M+\epsilon\, \mathrm{sin}M}$.

(1 )

Conversely, if the orbiting body is at aphelion at $\boldsymbol{t=0}$, then the Sun and the body are on opposite sides of the circle center at $\boldsymbol{t=0}$, and the shaded area in Fig. 3 must be added to the circular sector to find the area swept out by $\boldsymbol{\vec{r}(t)}$.  In this case,

$\large \boldsymbol{\theta (t)=M-\epsilon \, \mathrm{sin}M}$.

(2)

Equations 1 and 2 specify the position of the orbiting body as a function of time, relative to the center of the offset circle.  Both equations are needed below, as we will study the motion of two bodies that are either at perihelion or aphelion at $\boldsymbol{t=0}$.  Knowing $\boldsymbol{\theta (t)}$, it is straightforward to calculate the body’s position relative to the Sun.  Keep in mind that these results, and all of the calculations to follow, are accurate only to first order in $\boldsymbol{\epsilon }$.

### C.  Maintaining a Constant Separation between Spacecraft

Consider just the two “vertical” spacecraft A and B shown in Fig. 2.  Let the line $\boldsymbol{\overline{AB}}$ joining them be tilted by angle φ from the plane of the ecliptic, and suppose that the midpoint of $\boldsymbol{\overline{AB}}$ travels in a circular path about the Sun.  If $\boldsymbol{\phi < 90^{\circ}}$, then A and B are no longer equidistant from the Sun.  At the instant shown in Fig. 2, $\boldsymbol{r_{\mathrm{A,B}}=R\pm L\mathrm{cos}\phi /2}$.  To have a period of 1 yr, both spacecraft must move in elliptical orbits with semi-major axis $\boldsymbol{a=R=1\: \textrm{AU}}$.  If their orbital velocities are perpendicular to their position vectors at $\boldsymbol{t=0}$, A and B must be at aphelion and perihelion, respectively.  The aphelion and perihelion distances for an elliptical orbit are $\boldsymbol{r_{a,p}=a(1\pm \epsilon )}$, so $\boldsymbol{\epsilon =L\mathrm{cos}\phi /2R}$ for either orbit.  In addition, both orbits are inclined from the ecliptic plane by $\boldsymbol{\delta =\pm L\mathrm{sin}\phi /2R}$.  ($\boldsymbol{\left | \delta \right |< 0.5^{\circ}}$ for LISA).

Figure 4 shows the elliptical orbits of A and B as circles offset from the Sun by $\boldsymbol{\pm \epsilon R}$.  The two orbits still intersect, but now B passes through the intersection point ahead of A.  From Eqns. 1 and 2, when $\boldsymbol{M=\pi /2}$$\boldsymbol{\theta _{\mathrm{B}}-\pi /2=\pi /2-\theta _{\mathrm{A}}=\epsilon }$, so the separation s between the two spacecraft is

$\boldsymbol{s=2\epsilon R+R\, \mathrm{sin}(\theta _{\mathrm{B}}-\pi /2)+R\, \mathrm{sin}(\pi /2-\theta _{\mathrm{A}})\simeq 4\epsilon R=2L\mathrm{cos}\phi }$,

where we have used the small angle approximation $\boldsymbol{\mathrm{sin}\epsilon \simeq \epsilon }$ to simplify the expression.  Since the difference $\boldsymbol{\theta _{\mathrm{B}}-\theta _{\mathrm{A}}=2\epsilon\, \mathrm{sin}M}$, we expect $\boldsymbol{\left | s \right |}$ to oscillate between minima at $\boldsymbol{M=0,\pi }$ and maxima at $\boldsymbol{M=\pi /2,3\pi /2 }$ as the two spacecraft orbit the Sun.  This will be borne out in the following section.  But note that for $\boldsymbol{\phi =60^{\circ}}$, $\boldsymbol{s=L}$ at $\boldsymbol{M=0}$ and $\boldsymbol{M=\pi /2}$, suggesting that this may be the optimal angle of inclination for the plane of the space-based interferometer.

### D.  Numerical Calculation of Spacecraft Separation

Eqns. 1 and 2 allow us to calculate numerically the coordinates of spacecraft A and B, and find their separation, throughout their circumsolar orbits.  The calculations can be carried out using a common spreadsheet, such as Excel, which should make it easy for students to interpret and duplicate them.  First consider spacecraft A.  It is at aphelion at $\boldsymbol{t=0}$, and the plane of its orbit is tilted from the ecliptic plane by the small angle $\boldsymbol{\delta=L\, \mathrm{sin}\phi /2R}$.  Begin the spreadsheet by generating a column of values for the mean anomaly $\boldsymbol{0\leq M\leq 2\pi }$.  Next, use Eqn. 2 to calculate a corresponding column of values $\boldsymbol{\theta _{\mathrm{A}}=M-\epsilon\, \mathrm{sin}M}$, where $\boldsymbol{\epsilon =L\, \mathrm{cos}\phi /2R}$.  The coordinates of A, relative to the Sun, are then found using

$\large \boldsymbol{x_{\mathrm{A}}(t)=R(\mathrm{cos}\theta _{\mathrm{A}}+\epsilon )\mathrm{cos}\delta }$

$\large \boldsymbol{y_{\mathrm{A}}(t)=R\, \mathrm{sin}\theta _{\mathrm{A}}}$

$\large \boldsymbol{z_{\mathrm{A}}(t)=R(\mathrm{cos}\theta _{\mathrm{A}}+\epsilon )\mathrm{sin}\delta }$.

(3)

Spacecraft B is at perihelion at $\boldsymbol{t=0}$, so $\boldsymbol{\theta _{\mathrm{B}}(t)=M+\epsilon \, \mathrm{sin}M}$.  Its orbit in inclined by $\boldsymbol{-\delta }$, so its coordinates are given by

$\large \boldsymbol{x_{\mathrm{B}}(t)=R(\mathrm{cos\mathrm{\theta _{\mathrm{B}}}}-\epsilon )\mathrm{cos}(-\delta )}$

$\large \boldsymbol{y_{\mathrm{B}}(t)=R\, \mathrm{sin}\theta _{\mathrm{B}}}$

$\large \boldsymbol{z_{\mathrm{B}}(t)=R(\mathrm{cos\mathrm{\theta _{\mathrm{B}}}}-\epsilon )\mathrm{sin}(-\delta )}$.

(4)

After compiling the positions of A and B vs. time, their separation can be computed: $\boldsymbol{s(t)=\sqrt{(x_{\mathrm{A}}-x_{\mathrm{B}})^{2}+(y_{\mathrm{A}}-y_{\mathrm{B}})^{2}+\cdots }}$, and is shown in Figure 5 for tilt angles $\boldsymbol{\phi =0^{\circ},15^{\circ},\cdots ,90^{\circ}}$.  As anticipated, $\boldsymbol{s(t)=L}$ independent of time for $\boldsymbol{\phi =60^{\circ}}$.

For additional insight, it is useful to switch to a primed coordinate system rotating uniformly with angular speed $\boldsymbol{\Omega =2\pi /T}$ (so $\boldsymbol{\Omega t=M}$).  The $\boldsymbol{{x}'}$ axis points from the Sun through the midpoint of line $\boldsymbol{\overline{AB}}$, the $\boldsymbol{{y}'}$ axis points in the direction of increasing $\boldsymbol{\theta }$, and the $\boldsymbol{{z}'}$ axis is perpendicular to the plane of the ecliptic.  Transforming to the rotating coordinate system [Ref. 5],

$\large \boldsymbol{{x}'=x\, \mathrm{cos}M+y\, \mathrm{sin}M}$

$\large \boldsymbol{{y}'=-x\, \mathrm{sin}M+y\, \mathrm{cos}M}$

$\large \boldsymbol{{z}'=z}$.

(5)

Furthermore, by subtracting $\boldsymbol{R}$ from $\boldsymbol{{x}'}$, we can place the origin of the rotating system at the midpoint of line $\boldsymbol{\overline{AB}}$, which follows a circular path on the ecliptic.  Figure 6 shows the result for spacecraft A for $\boldsymbol{\phi =60^{\circ}}$.  Note that $\boldsymbol{{x}'-R}$,  $\boldsymbol{{y}'}$ and $\boldsymbol{{z}'}$ vary sinusoidally with period 1 yr, and that the value of $\boldsymbol{\sqrt{({x}'-R)^{2}+{y}'^{2}+{z}'^{2}}}$ is constant and equal to $\boldsymbol{L/2}$.  This indicates that as A orbits the Sun, it also executes a circular “epicycle” of radius $\boldsymbol{L/2}$ about the midpoint of $\boldsymbol{\overline{AB}}$, moving clockwise in the plane inclined by φ from the ecliptic.  Spacecraft B moves similarly, but occupies a position diametrically opposed to A on the same epicycle, so the separation between A and B is equal to $\boldsymbol{L}$ at all times.  If $\boldsymbol{\phi \neq 60^{\circ}}$, A and B still execute one epicycle per year, but it is elliptical rather than circular, with semi-axes $\boldsymbol{L/2}$ and $\boldsymbol{L\, \mathrm{cos}\phi }$.  In this case, the separation of the two spacecraft varies with time.

### E.  Maintaining the Shape of the Constellation

Imagine a fleet of n spacecraft arrayed about a central point that moves in a circular path about the Sun.  How should we choose their orbits so that their spacing and configuration are stable?  Generalizing from the preceding discussion, we want all of the spacecraft to move synchronously about a circular epicycle when viewed in the rotating coordinate system.  This is guaranteed if the fleet is embedded in a plane tilted from the ecliptic plane by $\boldsymbol{\phi =60^{\circ}}$.  If the spacecraft occupy the vertices of a regular polygon (equilateral triangle, square, pentagon, etc.), then each will pass through aphelion (or perihelion) at intervals $\boldsymbol{\Delta t=T/n}$ or $\boldsymbol{\Delta M=2\pi /n}$.  For example, in the preceding example of two spacecraft, B and A arrive at aphelion six months apart ($\boldsymbol{\Delta M=\pi }$).  Equally important, B’s orbit is the same as A’s, but rotated by 180° about the z-axis (perpendicular to the ecliptic plane).  So for n spacecraft, all orbits should have exactly the same shape and be inclined from the ecliptic by the same angle δ, but be rotated with respect to one another by $\boldsymbol{2\pi /n}$.

To illustrate, let us return to the four spacecraft shown in Fig. 2, and assume that $\boldsymbol{\phi =60^{\circ}}$.  Rather than following the ecliptic (as assumed in Section A), the orbits of C and D should have the same shape and tilt as A’s, but should be rotated about the z-axis by $\boldsymbol{\pm 90^{\circ}}$.  This is shown in Figure 7.  If A passes through aphelion at time $\boldsymbol{M=0}$, then C will pass through its aphelion at $\boldsymbol{M=\pi /2}$, while D has already done so at $\boldsymbol{M=-\pi /2}$.  At the instant shown in Fig. 2 ($\boldsymbol{M=0}$), C is moving upward through the ecliptic plane toward its aphelion, whereas D is moving downward.  Let’s check this by calculating the x and y positions of C and D at $\boldsymbol{M=0}$, which should be $\boldsymbol{x=R}$$\boldsymbol{y=\pm L/2}$.  From Eqn. 2, during the time interval $\boldsymbol{0< M< \pi /2}$, each spacecraft moves through an angle $\boldsymbol{\Delta \theta =\pi /2-\epsilon }$, measured from the center of its (offset) circular orbit.  Ignoring the small tilt of the orbits, the x and y coordinates (see Fig. 7) of C and D at $\boldsymbol{M=0}$ are, to first order in $\boldsymbol{\epsilon }$$\boldsymbol{x=R\, \mathrm{cos}\epsilon \simeq R }$ and $\boldsymbol{y=\pm (\epsilon R+R\,\mathrm{sin}\epsilon )\simeq \pm 2R\epsilon=\pm L/2 }$, as expected. [Ref. 7]

We conclude that to maintain its size and triangular configuration, LISA’s three spacecraft must lie in a plane tilted by 60° from the ecliptic plane, and must follow identical orbits rotated by 120°.  In fact, any number of spacecraft lying in this plane in any configuration will move synchronously if their orbits have the same semi-major axis a, eccentricity ε, and inclination δ, and in addition, are suitable rotated with respect to one another.  One should remember, though, that these results are valid only to first order in ε (or $\boldsymbol{L/R}$).  Exact solution of the equations of motion for LISA shows that the spacecraft separation will still vary by ± 1%. [Ref. 7]

### F.  Questions for Students

1. Draw a figure similar to Fig. 3, and derive Eqn. 2 for the position of a body that passes through aphelion at $\boldsymbol{t=0}$.

2.  A body is in an elliptical orbit about the Sun with period 1 yr and eccentricity  $\boldsymbol{\epsilon =0.10}$.  (a) What is the approximate eccentric anomaly θ two months after perihelion?  (b) Find the true anomaly ν and the distance r to the Sun at this time.  (Ans: (a) $\boldsymbol{M=60^{\circ}}$, so $\boldsymbol{\theta =65^{\circ}}$ (b) $\boldsymbol{\nu =70.4^{\circ}}$$\boldsymbol{r=0.96\: \textrm{AU}}$)

3. Use Fig. 1 to show that the eccentricity of LISA’s orbits is $\boldsymbol{\epsilon =L/(2R\sqrt{3})\simeq 0.005}$.

4. How can an L-shaped (as opposed to an equilateral triangle) fleet of three spacecraft be placed in an Earth-trailing orbit such that the spacing and configuration of the spacecraft remain stable?  (Ans: $\boldsymbol{a=1\: \textrm{AU}}$ for all three spacecraft,  which lie in a 60° plane; let $\boldsymbol{\epsilon =0}$ for the center craft, which follows the ecliptic around the Sun; $\boldsymbol{\epsilon =L/4R}$  for the other two, whose orbits are tilted by $\boldsymbol{\delta =\sqrt{3}L/4R}$ and are rotated by 90° with respect to each other.)

5. A fleet of five spacecraft is to be inserted in a circular orbit between Mars and Jupiter, 2.5 AU from the Sun.  The spacecraft are to remain in a regular pentagonal formation with minimal use of on-board rocket fuel.  (a) Explain how this can be done.  (b) One of the spacecraft passes through perihelion at $\boldsymbol{t=0}$.  How long before another passes through perihelion?  (Ans: $\boldsymbol{T=4\: \textrm{yr}}$, so $\boldsymbol{\Delta t=T/5=9.6\: \textrm{mo}}$.)

6. Suppose the plane of the fleet were inclined from the ecliptic by $\boldsymbol{\phi =120^{\circ}}$ (e.g., in Fig. 2, A is now closer to the Sun than B).  Would the configuration of the spacecraft be stable?

### Acknowledgments

We thank Jonathan Levine and Francesco Meschia for suggesting numerous ways to improve the clarity of the physics in this post.  The author is responsible for any remaining errors.

### References

1.  B. P. Abbott et al, “Observation of Gravity Waves from a Binary Black Hole Merger,”  Phys. Rev. Lett. 116, 061102 (2016)

2.  “LISA: Laser Interferometer Space Antenna: A cornerstone Mission for the Observation of Gravitational Waves,” ESA-SCI (2000) 11 July, 2000

3.  A nice treatment of a similar problem, by similar methods, is given by E. I. Butikov, “Relative Motion of Orbiting Bodies,” Am. J. Phys. 69, 63-67 (2001)

4.  This approximation was first used by ancient Greek astronomers to describe the motion of the Moon and Sun within the Ptolemaic (geocentric) system of the world.  See A. Berry, “A Short History of Astronomy,” (1898), available from Dover Publications.

5.  For an easy derivation of these transformation equations, see J. C. Amato and E. J. Galvez, Physics from Planet Earth, Chapter 2, problem 2.18 (Taylor & Francis/CRC Press, 2015).

6.  This topic is treated on a higher mathematical level by S. V. Dhurandhar et al, “Fundamentals of the LISA Stable Flight Formation,” Class. Quantum Grav. 22, 481-487 (2005).

7.  LISA Phase 0 CDF Study – Internal Final Presentation, ESTEC (European Space Research and Technology Centre) Conference, Noordwijk, Netherlands, March 8 – May 5, 2017,  slide 7, p. 110  http://sci.esa.int/lisa/59336-lisa-internal-phase-0-cdf-study-final-presentation/#

# The Range of North Korean ICBMs: Update

### A.  Introduction

This is an update of an earlier post (Post #13) published in August 2017.  On November 28, 2017, North Korea conducted a test flight of a heavier, longer-range missile (Hwasong-15) that appears to be capable of delivering a nuclear weapon anywhere in North America and Europe.  Moreover, additional information about all of North Korea’s missile tests is now publicly available, allowing a simpler and better estimate of missile range to be calculated.  In this post, we outline this new approach, and include an analysis of the Hwasong-15 test flight.

### B.  The “Effective Launch Speed”

The flight of an ICBM has three successive parts: the boost phase, while the missile’s rocket engines are firing; the unpowered midcourse phase, during which the missile is coasting under the influence of gravity; the terminal, or re-entry phase, when the warhead re-enters Earth’s atmosphere and converges on its target.  The midcourse phase occupies most of the flight time and is the easiest to observe.  Since August, the apogee (highest altitude) $\boldsymbol{h_{max}}$ reached in midcourse has been made public for each test flight.  With this new information, we can calculate an effective launch velocity and from it estimate the maximum range of the missile.  The numerical calculation of launch velocity based on total flight time, which we used in the previous post, is no longer necessary.

Recall that, for political reasons, North Korea uses highly “lofted” – nearly vertical  trajectories to test its ICBMs.  In the following, we will ignore the small horizontal component of the missile’s velocity when analyzing its test flight.  As before, we will also ignore the Earth’s rotation and neglect atmospheric drag.  (Air drag has little effect on the trajectory during the boost and midcourse phases because air density decreases rapidly with altitude, and throughout its boost phase, the ICBM is moving relatively slowly as it passes through the atmosphere.)

A missile’s midcourse trajectory is pre-determined by its burnout velocity $\boldsymbol{\vec{v}_{bo}}$ and altitude $\boldsymbol{h_{bo}}$ reached at the end of its boost phase.  Ignoring the small horizontal velocity of the missile, energy conservation requires

$\large \boldsymbol{\frac{E}{m}=\frac{1}{2}v_{bo}^{2}-\frac{Gm_{E}}{r_{E}+h_{bo}}=-\frac{Gm_{E}}{r_{E}+h_{max}}}$

(1)

where E is the total mechanical energy of the missile after burnout, m is the missile’s mass at burnout, G is the gravitational constant, and $\boldsymbol{m_{E}}$ and $\boldsymbol{r_{E}}$ are the mass and radius of Earth.  In the aftermath of a missile test, the only information publicly available are the total flight time and apogee $\boldsymbol{h_{max}}$: we do not know $\boldsymbol{v_{bo}}$ or $\boldsymbol{h_{bo}}$.  So imagine instead that the missile was launched at sea level by an instantaneous rocket burst (like an artillery shell), and it reached the same apogee as in the real flight test.  The “effective launch speed”  $\boldsymbol{v_{0}}$ due to this instantaneous burst obeys a relation akin to Eqn. 1,

$\large \boldsymbol{\frac{E}{m}=\frac{1}{2}v_{0}^{2}-\frac{Gm_{E}}{r_{E}}=-\frac{Gm_{E}}{r_{E}+h_{max}}}$.

(2)

For example, the Hwasong-14 missile launched on July 4 rose to a height of 2802 km (Ref. 1). Using Eqn. 2, its effective launch speed was $\boldsymbol{v_{0}=6.17\: \mathrm{km/s = 0.553 \boldsymbol{\, v_{esc}}}}$, where $\boldsymbol{v_{esc}=\sqrt{2Gm_{E}/r_{E}}}$ is the escape velocity (11.17 km/s) from the surface of the planet.

### C.  An Initial Estimate of Maximum Range

The previous section dealt only with near-vertical test flights.  Using those results, we can now estimate the range of the same missile when it is launched on a minimum energy trajectory (MET), i.e., the path that maximizes its range.  The physics is the same as we presented in Section C of the August post.  The missile trajectory is a segment of an ellipse with one focus ($\boldsymbol{F_{1}}$) located at the Earth’s center (Figure 2), and a semi-major axis a derived from conservation of energy.  In terms of the effective launch speed $\boldsymbol{v_{0}}$,

$\large \boldsymbol{\frac{E}{m}=\frac{1}{2}v_{0}^{2}-\frac{Gm_{E}}{r_{E}}=-\frac{Gm_{E}}{2a}}$.

(3)

Solving for a and using $\boldsymbol{u_{0}=v_{0}/v_{esc}}$, we obtain

$\large \boldsymbol{\frac{2a}{r_{E}}=\frac{1}{1-u{_{0}^{2}}}}$.

(4)

In Figure 2, the launch point is P and the target point is Q, $\boldsymbol{F_{2}}$ is the second (empty) focus, and θ is the angle between $\boldsymbol{\overline{F_{1}P}}$ (or $\boldsymbol{\overline{F_{1}Q}}$) and the major axis of the ellipse.  The missile’s range is maximized when θ is a maximum, which occurs when $\boldsymbol{\overline{F_{2}P}}$ is perpendicular to the major axis.  Since $\boldsymbol{\overline{F_{1}P}+\overline{F_{2}P}=2a}$,

$\large \boldsymbol{sin(\theta _{max})=\frac{2a-r_{E}}{r_{E}}=\frac{u_{0}^{2}}{1-u_{0}^{2}}}$

(5)

and the maximum range is

$\large \boldsymbol{R_{max}=2r_{E}\theta _{max}}$.

(6)

To illustrate, for the Hwasong-14 launched on July 4, we found that $\boldsymbol{u_{0}=0.553}$, so $\boldsymbol{\theta _{max}=0.455\: \mathrm{(rad)}=26.1^{\circ}}$ and  $\boldsymbol{R_{max}=5810\, \mathrm{km}}$.  However, this estimate is significantly shorter than the published range of 6700 km (Ref. 2), and we will correct it in Section D.

Using the reflection property of the ellipse, the launch angle (measured relative to the local horizontal), is given by

$\large \boldsymbol{\varphi _{0}=\frac{\pi }{4}-\frac{\theta _{max}}{2}}$.

(7)

To illustrate, for the July 4 missile test, $\boldsymbol{\theta _{max}=0.455\: \mathrm{(rad)}}$, so $\boldsymbol{\varphi _{0}=.558\: \mathrm{(rad)}=32.0^{\circ}}$.  See the August post for more details.

### D.  An Improved Estimate of Missile Range

The range calculated using the above procedure  will be consistently shorter than the figure quoted by arms control experts, because the launch velocity  calculated for the near-vertical test flight is not equal to $\boldsymbol{v_{0}}$ for the same missile on a MET.  Let’s look at this in more detail.  In a typical ICBM flight, the missile first rises vertically from its launch pad (or silo), then is tipped by a few degrees toward its target by firing auxiliary steering thrusters or, in more sophisticated designs, by briefly diverting (gimbaling) the main engine exhaust off-axis.  The rocket thrust is then re-aligned with the missile axis, which then rotates to the desired launch angle $\boldsymbol{\varphi _{0}}$ using gravity to alter its direction.  This so-called “gravity turn” minimizes stress on the missile body, and is a good approximation to what actually occurs in practice (Ref. 3).  Figure 3 shows boost phase trajectories of a missile undergoing a gravity turn for several values of $\boldsymbol{\varphi_{0}}$, plus a free body diagram for the missile .  The equations of motion for the missile during this maneuver are:

$\large \boldsymbol{a_{\parallel }=\frac{dv}{dt}=\frac{F_{th}}{m}-gsin\varphi }$

$\large \boldsymbol{a_{\perp }=v\frac{d\varphi }{dt}=-gcos\varphi }$

(8)

where $\boldsymbol{F_{th}}$ is the rocket thrust, m is the instantaneous missile mass, and φ is the instantaneous angle  ($\boldsymbol{0\leq \varphi \leq \varphi _{0}}$) between the velocity vector $\boldsymbol{\vec{v}}$ and the horizontal direction.  $\boldsymbol{a_{\parallel }}$ and $\boldsymbol{a_{\perp }}$ are the components of the acceleration parallel and perpendicular to $\boldsymbol{\vec{v}}$.  The “gravity loss” term $\boldsymbol{-gsin\varphi }$ reduces the burnout velocity, and its integrated effect increases with the launch angle $\boldsymbol{\varphi _{0} }$ and the duration of the boost phase. Since $\boldsymbol{F_{th}}$, mφ, and $\boldsymbol{g=g(h)}$ are functions of time, the two (coupled) Eqns. 8 above must be solved numerically to find the burnout speed $\boldsymbol{v_{bo}}$ and altitude $\boldsymbol{h_{bo}}$  (Fig. 3).  Clearly, the more we know about a missile’s design, the more accurately we can model its boost phase behavior.

We want to find the $\boldsymbol{\varphi _{0}}$-dependence of $\boldsymbol{v_{bo}}$ and $\boldsymbol{h_{bo}}$, and thereby obtain a $\boldsymbol{\varphi _{0}}$-dependent correction to the effective launch speed $\boldsymbol{v_{0}}$.  While the design details of North Korea’s Hwasong-14 and Hwasong-15 ICBMs are secret, they are known to be two-stage liquid-fueled rockets.  So let us instead study the boost phase behavior of two similar missiles of known design.  One well-characterized missile is the Titan-II, a silo-based American ICBM of 1960s vintage (Ref. 4).  A second missile design comes from a study of the feasibility of ballistic missile defense, conducted by the American Physical Society in 2004 (Ref. 5).  This missile, called Model L in the study, has a significantly higher thrust and shorter burnout time than the Titan-II, resulting in a higher burnout velocity.  Remarkably, in spite of their different designs, the $\boldsymbol{\varphi _{0}}$-dependent correction to $\boldsymbol{v _{0}}$ is almost exactly the same for both missiles.

Using Eqns. 8, we simulated the boost phase of both missiles for seven values of $\boldsymbol{\varphi _{0}}$  .  In each simulation, the missile left its launch pad and traveled straight upwards for 10 s before being given a small angular “kick” to begin the gravity turn.  The kick angle (< 5°) was selected by trial and error to produce the desired launch angle  $\boldsymbol{\varphi _{0}}$ at burnout ($\boldsymbol{11.25^{\circ}\leq \varphi _{0}\leq 90^{\circ}}$).  For each missile and launch angle, $\boldsymbol{v_{bo}}$ and $\boldsymbol{h_{bo}}$ were recorded, and an equivalent launch speed $\boldsymbol{v_{0}}$ was calculated by combining Eqns. 1 and 2:

$\large \boldsymbol{v_{0}^{2}=v_{bo}^{2}+2Gm_{E}\left ( \frac{1}{r_{E}}-\frac{1}{r_{E}+h_{bo}} \right )}$.

(9)

For example, in Section C we found $\boldsymbol{v_{0}(90^{\circ})=6.17\: \mathrm{km/s}}$ and $\boldsymbol{\varphi _{0}=32^{\circ}}$ in our analysis of of the July 4 missile test.  According to Fig. 4, we should add $\boldsymbol{\Delta v_{0}=0.31\pm 0.01\: \mathrm{km/s}}$ to $\boldsymbol{v_{0}(90^{\circ})}$ to improve our estimate of range.  Using Eqn. 5, the corrected launch speed $\boldsymbol{{v}'_{0}=6.48\: \mathrm{km/s}}$ yields a new value $\boldsymbol{{\theta }'_{max}=30.9^{\circ}}$ and a new range estimate $\boldsymbol{{R}'_{max}=6790\pm 40\: \mathrm{km}}$.  This agrees well with the estimate (≈ 6700 km) adopted by arms control analysts.

The new value $\boldsymbol{{\theta }'_{max}}$ changes the value of $\boldsymbol{\varphi _{0}}$, which in turn influences $\boldsymbol{\Delta v_{0}}$.  This suggests that we iterate our calculations until they converge.  We think this is unwarranted, however, because we have introduced an unavoidable error in $\boldsymbol{\Delta v_{0}}$ by replacing the Hwasong missiles with the Model L or Titan-II.

Table 1 summarizes our results for the North Korean ICBM tests of 2017, and compares them against the range estimates appearing in the news media.  Also included is the May 14 test of a Hwasong-12 intermediate range ballistic missile (IRBM).  Note that the uncertainty increases with launch speed: as $\boldsymbol{v_{0}\rightarrow 7.90\: \mathrm{km/s}}$, which is the speed needed for low Earth orbit (i.e., unlimited range), our range estimate becomes more sensitive to errors in $\boldsymbol{v_{0}}$.  Perhaps, for a similar reason, the maximum range of the Hwasong-15 has only been published as a lower limit.

### References

1. Ankit Panda, “The Hwasong-15: Anatomy of North Korea’s New ICBM,” The Diplomat, December 6, 2017, https://thediplomat.com/2017/12/the-hwasong-15-the-anatomy-of-north-koreas-new-icbm/

2. Published estimates of missile range are originally due to D. Wright, All Things Nuclearhttp://allthingsnuclear.org/dwright/ , May 13, July 3, July 28, September 14 and November 28, 2017

3.  Glen J. Culler and Burton D. Fried, “Universal Gravity Turn Trajectories,” J. Appl. Phys. 28, 672-676 (1956)

4.  D. K. Stumpf, “Titan-II: A History of a Cold War Missile Program,” U. Arkansas Press (2000), Chapter 8

5.  D. K. Barton et al, “Report of the APS Study Group on Boost-Phase Intercept Systems for National Missile Defense,” Rev. Mod. Phys76, S1 (2004)

# Studying Haumea: A Dwarf Planet with Two Moons and a Ring

## A.  Introduction

Haumea is a Pluto-sized dwarf planet lying beyond Neptune within the Kuiper Belt (Fig. 1).  Discovered in 2004, it is the fifth largest dwarf planet known, and is unique among them because it hosts a ring in addition to its two moons.  It is also one of the fastest spinning objects in the solar system, completing one revolution in less than 4 hours.  The ring was discovered in 2017 by Ortiz et al (Ref. 1) when Haumea passed between Earth and a distant star, temporarily blocking its light from reaching Earth.  This stellar occultation was detected at nine observatories in central and eastern Europe, and the light curves recorded at these sites were used to determine the ring’s radius and opacity in addition to the dimensions of the dwarf planet itself.  The Haumea system offers a wealth of opportunity for illustrating basic physics concepts to first year college students.  Lying in the outer reaches of the solar system, its physical properties are likely to be unknown a priori (in contrast to Saturn, for example), making it an ideal topic for study or testing.  In the following, we offer a simple treatment of the stellar occultation data, plus an assortment of problems suitable for homework exercises, quizzes or exams in an introductory physics course.  We hope you will find these materials useful in your teaching.

## B.  The Size and Shape of Haumea

The occultation occurred on January 21, 2017.  Figure 2 shows the light curves recorded at the nine observatories participating in the study.  The large dip centered near 200 s is due to occultation by the dwarf planet, and the narrower, weaker dip near 50 s is due to the ring.  Figure 3 shows the relative motion of Haumea and the occulted star as observed at each observatory.  The blue lines are the occultation chords, i.e., the apparent motion of the star as the planet passed in front of it, as seen at each observatory.  Note that the length of the chord observed at Konkoly (in Hungary) is nearly equal to the diameter b-b′ of Haumea.

To use the light curves to determine the chord lengths, we need to know the proper speed of Haumea relative to Earth.  The proper speed is the component of the relative velocity perpendicular to the line of sight from Earth to the dwarf planet.  Using the online orbital calculator on JPL’s Horizons website (Ref. 2), it is straightforward (but tedious) to find the proper speed at the time of occultation: $v_{\perp }=14.367\textup{ km/s}$.

1.  The ingress/egress times recorded at Konkoly are 3:08:20.3 and 3:10:17.4 UTC.      Find the length of the Konkoly occultation chord (in km) shown in Fig. 3.  (Ans:  1680 km)

In the following questions, assume that the shape of Haumea is a triaxial ellipsoid with semi-axes a > b > c.

2.  From direct measurements on Fig. 3, the “diameter” b-b′ is 2% greater than the Konkoly chord.  Calculate the semi-axis b of the dwarf planet.  (Ans: 857 km.  Compare to the published value 852 ± 4 km)

3.  Haumea rotates rapidly around its shortest axis, labeled c in Fig. 3.  At the time of the occultation, the brightness of the planet observed from Earth was at its minimum value.  Its maximum brightness is 34% greater than this minimum.  Find the length of the semi-axis a.  (Ans: 1150 km.  Compare to the published value 1161 ± 30 km.)

4.  The length of the rotation axis c = 513 ± 16 km, found from analysis of the brightness vs. time curve.  The volume of a triaxial ellipsoid is $V=\frac{4}{3}\pi abc$.  Using the dimensions found above, calculate the radius R of a sphere having the same volume.  (Ans: 794 km.  The published value is 797.5 ± 5.5 km.)

5.  Phil Plait, who writes the Bad Astronomy blog (Ref. 3), likened the shape of Haumea to a “flattened potato.”  Why is the planet flattened?  (See PPE, section 8.14.)

## C.  The Ring and Moons of Haumea

The ring surrounding Haumea is in the plane of its equator.  According to the online Extended Data, Table 2 of Ref. 1, occultation by the ring was seen at Ondrejov (Czech Republic) at 3:06:48.4 UTC, while the center of the occultation by the dwarf planet occurred at 3:09:20.7 UTC.  (See Figs. 2 and 4.)

6.  Estimate the radius of the ring.  (Ans: $2188\textup{ km}$.  The published value is $2287_{+75}^{-45}\textup{ km}$.)

In the following questions, assume that Haumea is a spherically symmetric body with the radius R given in Question 4 above.

7.  Haumea’s largest moon, Ha’iaka, has a near-circular orbit in the equatorial plane of the planet (Ref. 4).  Its period is $P=49.462\pm .083\textup{ d}$, and its semi-major axis is $a=49,880\pm 200\textup{ km}$.  Find the mass of Haumea.  (Ans: $4.02\times 10^{21}\textup{ kg}$)

8.  A second moon, Namaka, orbits Haumea with period ${P}'=18.278\pm .008\textup{ d}$.  Find the semi-major axis of its orbit.  (Ans: ${a}'=25,680\textup{ km}$)

9.  Calculate the average density of Haumea, and compare it to that of Pluto (Ref. 5), another dwarf planet lying within the Kuiper Belt.  Pluto has a radius of 1151 km and a mass of $1.309\times 10^{22}\textup{ kg}$.  (Ans: $\rho _{H}=1.885\times 10^{3}\mathrm{\: kg/m^{3}}$, $\rho _{P}=2.049\times 10^{3}\mathrm{\: kg/m^{3}}$.  Nearly the same!)

10.  Show that Haumea’s moons lie well beyond its Roche limit. (See PPE, section 8.16.) (Ans: $R_{Roche}=1950\textup{ kg}$.)

Note:  The Roche limit derived in PPE, Eqn. 8.22, assumes that the density of the main body and its satellites are the same.  If Haumea has a rocky core and an icy crust, as suspected from its albedo, and if its ring and moons were formed by impact, then they are likely to have a density closer to that of ice, $\rho\approx 1000\mathrm{\: kg/m^{3}}$.  As shown in the discussion preceding Eqn. 8.22 in PPE, this lower density would increase the magnitude of the Roche limit.  Ortiz et al assume a satellite density of $500\mathrm{\: kg/m^{3}}$ to obtain a Roche limit substantially larger than the ring radius.

## D.  Surface Gravity on Haumea

In the following exercises, assume once again that Haumea is a spherically symmetric body of radius 797.5 km and density $1.885\times 10^{3}\mathrm{\: kg/m^{3}}$, and that its period of rotation is 3.9155 ± .0001 hr.

11.  What is the gravitational acceleration $g$ on the surface of Haumea?  What is the escape velocity from the surface?  (Ans: $0.42\mathrm{\: m/s^{2}}$, 820 m/s)

12.  Picture a tiny object lying on the surface of the spinning planet.  Draw a free body diagram for the object.  At what location(s) is the normal force acting on it a maximum?  Where is the normal force a minimum?  (Ans: maximum at the poles, minimum at the equator)

13.  If the dwarf planet is a gravitationally bound “rubble pile,” it will begin to break apart when the minimum normal force on a surface object approaches zero.  (See PPE, problem 8.15.)  What is the minimum period of rotation for Haumea to remain stable?  (Ans: 2.41 hr)

14.  A simple pendulum of length 1 m is placed on the surface of Haumea.  What is its period of oscillation if it is located (a) at the north pole; (b) on the equator?  (Ans: (a) 9.70 s  (b) 12.3 s)

15.  Imagine that the pendulum is located at latitude $\theta$ and is not swinging.  Upon close inspection, you find that it is not hanging vertically, but is deflected by an angle $\phi$.  Derive an equation showing the dependence of $\phi$ on the latitude and on Haumea’s angular velocity $\omega$.  Is the pendulum bob deflected toward the pole or toward the equator?  Justify this with a free body diagram.  What is the deflection angle when $\theta =45^{\circ}$ degrees?

(Ans: $\large \boldsymbol{tan\phi =\frac{\omega ^{2}Rsin2\theta }{2\left ( g-\omega ^{2}Rcos^{2}\theta \right )}}$, toward the equator, $\large \phi (45^{\circ})=13.1^{\circ}$)

## E. References

1.  J. L. Ortiz et al, Nature 550, 219 (October 12, 2017)

4.  D. Ragozzine and M. E. Brown, Astro. J. 137, 4766 (June 2009)

5.  See, e.g., https://solarsystem.nasa.gov/planets/pluto/facts

# The Last Days of Cassini: the Grand Finale

### A.  Introduction

A central premise of Physics from Planet Earth is that students will be more motivated to study physics if their coursework addresses important and newsworthy scientific events.  A good example is the Cassini-Huygens mission to Saturn, which ended dramatically on September 15, 2017 when the spacecraft dipped into the Saturnian atmosphere and self-destructed — its grand finale.  Cassini‘s fiery demise was carefully planned by NASA seven years ago (Ref. 1), when the spacecraft’s fuel supply began to run low, limiting opportunities for further exploration of Saturn and its surroundings.  Flight engineers chose to end the mission by crashing into Saturn to eliminate any possibility of the spacecraft contaminating Enceladus or Titan, two moons of Saturn that potentially harbor life.

Cassini was launched from Cape Canaveral 20 years ago.  After 4 gravity assists (2 from Venus, 1 from Earth, and 1 from Jupiter), it entered orbit about Saturn in July 2004, the first spacecraft ever to do so.  The ESA’s (European Space Agency) Huygens spaceprobe was dispatched from Cassini on Christmas day, 2004, and soft-landed on Titan three weeks later, the most distant landing ever achieved by a spacecraft from Earth.  Over the next 13 years, Cassini discovered 7 new moons of Saturn (bringing the total to 62) and conducted close flybys of 13 moons.  To visit them, it employed gravity assists, or “targeted flybys,” to alter its flight trajectory with minimal expenditure of fuel.  Most of these flybys — 127 in all — involved Titan, Saturn’s largest moon and the only one with sufficient mass to alter Cassini‘s trajectory significantly.  By taking advantage of gravity assists, NASA flight engineers were able to extend the mission by 9 years.  Their efforts were amply rewarded.  Together, Cassini and Huygens revolutionized our understanding of the Saturnian system.  In addition to providing stunning close-up images of the planet, its rings and moons, the two spacecraft discovered subsurface water oceans on Titan, Enceladus and Dione, surface oceans of methane (CH4) on Titan, and gaseous plumes  containing H2O, CO2, CH4 and H2 emanating from Enceladus, making the latter the most promising candidate in our solar system for hosting extraterrestrial life.

Prior to its “grand finale” orbits, Cassini used a targeted flyby of Titan (designated T125) to enter a “ring-grazing” orbit, with periapse just outside Saturn’s F ring.  After completing 20 ring-grazing orbits, the spacecraft executed its final targeted flyby (T126) of Titan on April 22, 2017, and was inserted into a daring trajectory that — near its periapse — passed between the cloud tops of Saturn and its innermost ring, a gap only 3000 km wide.  Over the next 5 months, Cassini executed 22 grand finale orbits, until a high altitude flyby of Titan on September 12 caused the spacecraft to plunge into Saturn’s atmosphere three days later.

In this post, we focus on T126 to learn how Cassini was transferred from its last ring-grazing orbit to the first of its grand finale orbits.  This is done via a set of homework exercises and/or exam problems suitable for first year university physics students.  We gratefully acknowledge the generous help of Dr. Duane Roth, Cassini Navigation Team Chief, who kindly forwarded Reference 1 and provided detailed instructions for using the online software package WebGeoCalc (Ref. 2) to calculate the trajectory of Cassini relative to Saturn and Titan, both before and after T126.  We drew heavily WebGeoCalc’s calculations to write this post.  Thank you, Dr. Roth!

Cassini‘s grand finale was a headline-grabbing event, and perfectly timed for inclusion in the fall semester syllabus of an introductory mechanics course.  We hope your students will enjoy learning the physics behind the headlines.  Farewell, Cassini!  Right up to your final moments, you performed brilliantly!

### B.  The T126 flyby

In the exercises below, the T126 flyby should be treated as an elastic “collision” between Cassini and Titan.  The following information will be useful for solving the problems.

Titan is in a near-circular orbit (ε = 0.0288) about Saturn, well-aligned to the planet’s equatorial plane, with period $\boldsymbol{P_{T}=1.379\times 10^{6}}$ s (16.0 d).  Using Kepler’s 3rd law, its semi-major axis is found to be $\boldsymbol{a_{T}=1.222\times 10^{6}}$ km.  Prior to the T126 flyby, Cassini was in a highly eccentric orbit (ε = 0.795) about Saturn with inclination about 60° relative to the equatorial plane, period $\boldsymbol{P_{C}=6.205\times 10^{5}}$ s (7.2 d) and semi-major axis $\boldsymbol{a_{C}=7.178\times 10^{5}}$ km.  The flyby occurred at a distance $\boldsymbol{r=1.206\times 10^{6}}$ km from Saturn, not far from the spacecraft’s apoapse.  (See Figure 2.)  At the encounter point, the spacecraft crossed the plane of Titan’s orbit from below to above.  As a result of the flyby, Cassini entered an even more eccentric orbit (ε = 0.906) which, at periapsis, passed between the planet and its innermost ring.

Exercises for students:

1. Titan is in a near-circular orbit about Saturn with period $\boldsymbol{P_{T}=1.379\times 10^{6}}$ s.  The T126 flyby occurred when Titan was $\boldsymbol{1.206\times 10^{6}}$ km from Saturn.  Find its speed at this time.  Use $\boldsymbol{GM_{S}=3.793\times 10^{7}\:\mathrm{km^{3}/s^{2}}}$.  (Ans: $\boldsymbol{v_{T}=5.645\: \mathrm{km/s}}$)

2.  Just before the flyby, Cassini was $\boldsymbol{1.19\times 10^{6}\: \mathrm{km}}$ from Saturn and moving with speed $\boldsymbol{v_{C}=3.291\: \mathrm{km/s}}$.  The eccentricity of its orbit was 0.795.  Show that, at periapsis, the spacecraft was just outside Saturn’s F ring, which has a radius of about $\boldsymbol{1.40\times 10^{5}\: \mathrm{km}}$.

3.  Titan’s orbit is within 0.3° of the plane of Saturn’s rings, whereas Cassini‘s orbit is tilted by about 60°.  As shown in this video,  Entering Final Orbits: Cassini’s Grand Finale (Ref 3), Cassini passed in front of Titan during T126.  (The encounter occurs 20 s from the beginning of the video.)  Let $\boldsymbol{\vec{u}}$ and $\boldsymbol{\vec{{u}'}}$ be the velocities of the spacecraft relative to Titan before and after the flyby.  Make a rough sketch of the velocity vectors $\boldsymbol{\vec{v}_{T}}$$\boldsymbol{\vec{u}}$ and $\boldsymbol{\vec{{u}'}}$, and explain why the spacecraft’s speed decreased due to the flyby.  No calculations are necessary.

4.  In a reference frame whose x– axis is parallel to $\boldsymbol{\vec{v}_{T}}$ , and whose z-axis is perpendicular to the moon’s orbit, the velocity of Titan just before the flyby may be expressed as $\boldsymbol{\vec{v}_{T}=5.645\, \hat{i}\: \mathrm{km/s}}$, and the velocity of Cassini as $\boldsymbol{\vec{v}_{C}=1.161\, \hat{i}-0.701\, \hat{j}+2.998\, \hat{k}\: \mathrm{km/s}}$.   Calculate the angle between the two vectors and also find the relative velocity vector $\boldsymbol{\vec{u}=\vec{v}_{C}-\vec{v}_{T}}$.  (Ans: 66°)

5.  As Cassini approached Titan, the angle between the relative velocity vector $\boldsymbol{\vec{u}}$ and $\boldsymbol{\vec{v}_{T}}$ was 145.5°.  The flyby deflected the direction of the relative velocity by 9.1°.  Show that this reduced the spacecraft’s speed relative to Saturn by about 850 m/s.

6.  True or false:  the change in the speed of the spacecraft is equal to $\boldsymbol{\left | {\vec{v}}'_{C}-\vec{v}_{C} \right |}$, i.e., the absolute value of the change in the velocity of the spacecraft.

7.  At its closest approach to Titan, Cassini was only 980 km above the surface of the moon, or 3554 km from its center.  What was its speed relative to Titan at this time?  (Ans:  $\boldsymbol{u_{max}= 5.880 \: \mathrm{km/s}}$)

8.  The unfueled mass of Cassini was 2125 kg, and the exhaust speed of its MMH/NTO fuel mixture is 2.8 km/s.  What is the minimum mass of fuel that would have been needed to change the spacecraft’s speed (relative to Saturn) by 850 m/s?  Assume that the maneuver totally depleted the fuel supply.  (Ans: 750 kg)

9.  In its new trajectory (after T126), Cassini made 22 full orbits with period $\boldsymbol{5.565\times 10^{5}\: \mathrm{s}}$ before a long distance encounter with Titan caused it to plunge into the Saturnian atmosphere.  Why 22 orbits?  (Hint: compare the periods of Cassini and Titan.)

The following two exercises can be completed after studying Chapter 12 of PPE or similar material.

10.  Prove that in the reference frame of Titan, the spacecraft was deflected by 9.1° due to T126.  (This is the angle between $\boldsymbol{{\vec{u}}'}$ and $\boldsymbol{\vec{u}}$.

11.  On September 12, Cassini came within $\boldsymbol{1.2\times 10^{5}\: \mathrm{km}}$ of Titan, and its relative speed at closest approach was 5.2 km/s.  Estimate the angle by which the spacecraft was deflected (in Titan’s reference frame), and the maximum change of speed relative to Saturn.  This change in speed was enough to ensure that the spacecraft entered Saturn’s atmosphere at periapse.

### References

1.  Brent Buffington et al, IAC-10-C1.9.2, 61st International Astronomical Congress, Prague, CZ (2010)

3.  Entering Final Orbits: Cassini‘s Grand Finale.  Video is available courtesy of NASA/JPL at https://saturn.jpl.nasa.gov/resources/7508/?category=videos .

# The Range of North Korean Ballistic Missiles

### A. Introduction

North Korea’s missile launches of July 2017  demonstrate that it has developed ballistic missiles (Figure 1) with intercontinental range (ICBMs).  This, plus its parallel development of nuclear weapons, has added a grave new threat to world order.  This is a grim topic, but because it has such immediate importance — and encompasses fundamental physical principles — it deserves to be discussed in our introductory physics classes.

The missile launch of July 4 flew for 37 minutes (Ref. 1), and the one launched on July 28 was aloft for 47 minutes (Ref. 2).  Based on this information, analysts in the United States and Japan estimate the maximum range of the former to be 6,700 km, and the latter to be 10,400 km, placing Los Angeles, Denver, Chicago, New York, Boston and other North American population centers within targeting range of Pyongyang.  The purpose of this post is to show, using mathematics and concepts appropriate to the introductory calculus-based mechanics course, how to estimate missile range from the (observable) flight time.  Following a treatment of Kepler’s Laws, it should be possible to present this material in a single class period.

Note:  See also Post 15: The Range of North Korean Ballistic Missiles: Update (Jan 14, 2018) for a better approach to this issue.  In particular, the numerical technique for estimating the burnout speed from the total flight time is no longer necessary.

### B.  Determining the Burnout Speed of a Ballistic Missile

To simplify the discussion, assume the missile altitude at burnout (when rocket fuel is exhausted) is negligible compared to the Earth’s radius, neglect atmospheric drag, and also neglect the Earth’s rotation.  North Korea uses highly “lofted” trajectories to test its ICBMs; that is, they are launched at much steeper angles (approaching 90º) than they would be to achieve maximum range.  In this section, we will assume that the missile is launched vertically, and calculate its flight time $\boldsymbol{t_{f}}$ and maximum altitude $\boldsymbol{h_{max}}$ as functions of its burnout speed $\boldsymbol{v_{0}}$.  For small $\boldsymbol{v_{0}}$ and $\boldsymbol{h_{max}}$, this is an easy task, as the gravitational acceleration $\boldsymbol{g}$ is roughly constant, and the familiar equations of kinematics are applicable: $\boldsymbol{t_{f}=2v_{0}/g}$ and $\boldsymbol{h_{max}=\frac{1}{8}gt_{f}^{2}}$.  (Note: $\boldsymbol{t_{f}}$ is twice the time from launch to maximum altitude.)  But a lofted ICBM will reach an altitude comparable to the Earth’s radius $\boldsymbol{r_{E}}$, so the gravitational acceleration varies significantly during the missile’s flight.  This complicates the analysis considerably.

Let $\boldsymbol{G}$ be the gravitational constant, $\boldsymbol{M}$ be the mass of the Earth, and $\boldsymbol{r}$ be the instantaneous distance of the missile from the center of the planet.  Launched vertically, the speed of the missile at maximum altitude ($\boldsymbol{r=r_{max}}$) is zero, so by energy conservation,

$\large \boldsymbol{\frac{1}{2}v_{0}^{2}-\frac{GM}{r_{E}}=0-\frac{GM}{r_{max}}}$

or

$\large \boldsymbol{r_{max}=\frac{r_{E}}{1-r_{E}v_{0}^{2}/2GM}=\frac{r_{E}}{1-u_{0}^{2}}}$

(1)

where $\boldsymbol{u_{0}=v_{0}/v_{esc}}$ and $\boldsymbol{v_{esc}=\sqrt{2GM/r_{E}}}$ is the escape velocity from the planet’s surface.

Our first task is to find the flight time $\boldsymbol{t_{f}}$ in terms of the burnout speed $\boldsymbol{u_{0}}$.  One way to proceed is to integrate the equation of motion numerically.  The most straightforward numerical technique is the second-order Taylor approximation (Ref. 3)in which the stepwise equations for position and speed are reminiscent of the kinematic equations for projectile motion:

$\large \boldsymbol{r_{n+1}=r_{n}+v_{n}\Delta t+\frac{1}{2}a_{n}\Delta t^{2}}$

and

$\large \boldsymbol{v_{n+1}=v_{n}+\frac{1}{2}\left ( a_{n}+a_{n+1} \right )\Delta t}$

where the acceleration $\boldsymbol{a_{n}=-GM/r_{n}^{2}}$ and $\boldsymbol{\Delta t}$ is the time step.  A simple Excel spreadsheet is adequate for the integration.  Figure 2 shows the computed altitude $\boldsymbol{h(t)=r(t)-r_{E}}$ for a projectile that has reached its maximum altitude at $\boldsymbol{t=0}$, for the three values $\boldsymbol{u_{0}}$ = 0.5, 0.6 and 0.7.  The intersection of each curve with the x-axis ($\boldsymbol{h=0}$) marks half the total flight time of the missile.

Some instructors may prefer an analytic expression for $\boldsymbol{h(t)}$.  A power series solution for $\boldsymbol{h(t)}$, due to Foong (Ref. 4), is outlined in the Appendix and is displayed in Fig. 2 along with the numerical solution.  Figure 3 shows the total flight time $\boldsymbol{t_{f}}$ as a function of $\boldsymbol{u_{0}}$, calculated numerically as well as by the power series approximation.  Knowing $\boldsymbol{t_{f}}$, one may find $\boldsymbol{u_{0}}$ directly from Fig. 3.

### C.  Finding the Maximum Range of the Missile

When the burnout speed is low ($\boldsymbol{u_{0}\ll 1}$), the maximum altitude and range of the missile are small, and its trajectory is adequately described as a parabola over a flat Earth.  At the higher speeds needed for intercontinental targeting, the parabolic approximation is no longer adequate, and the trajectory must be treated properly as a segment of an ellipse having one focal point at the center of the planet.  (This follows from Kepler’s Laws.) Changing the launch angle changes the eccentricity $\boldsymbol{\varepsilon }$ of the ellipse  — which affects the range of the missile — but not the length of its semi-major axis $\boldsymbol{a}$.  The goal of this section is to find the optimum launch angle (the one that yields the maximum range), and to calculate that range as a function of the burnout speed $\boldsymbol{u_{0}}$.  This can be done easily by appealing to the reflection property of an ellipse (Ref. 5).

In Figure 4, line $\boldsymbol{\overline{mn}}$ is tangent to the ellipse at point P.  The reflection, or focal property states that $\boldsymbol{\overline{F_{1}P}}$ and $\boldsymbol{\overline{F_{2}P}}$ make equal angles with the tangent line.  So at point P, angle φ = φ′ , and at point P′, angle γ = γ′.  If the ellipse were a reflective surface, any light ray emanating from one focus would, after reflection, pass through the other focus.  For example, a parabola is an ellipse with one focus at infinity.  As is well known, light rays originating at infinity (parallel rays) converge at the focal point of a parabolic mirror.

Next we review the energetics of orbital motion governed by an inverse square central force.  If the eccentricity ε of the elliptical orbit is zero, the trajectory is circular with radius r.  The orbital speed is found by equating the centripetal acceleration to the gravitational force per unit mass, $\boldsymbol{v^{2}/r = GM/r^{2}}$, so

$\large \boldsymbol{\frac{E}{m}=\frac{1}{2}v^{2}-\frac{GM}{r}=-\frac{GM}{2r}}$

where is the total mechanical energy of the missile and m is its mass.  When $\boldsymbol{\varepsilon \neq 0}$, a similar relation holds (Ref. 6): the energy per unit mass is fixed — not by the orbital radius r but by the length of the semi-major axis a:

$\large \boldsymbol{\frac{E}{m}=\frac{1}{2}v^{2}-\frac{GM}{r}=-\frac{GM}{2a}}$.

An ICBM’s energy and semi-major axis are set by its burnout speed $\boldsymbol{v_{0}}$, which we have assumed is reached when the missile is close to the Earth’s surface ($\boldsymbol{r\simeq r_{E}}$):

$\large \boldsymbol{\frac{E}{m}=\frac{1}{2}v_{0}^{2}-\frac{GM}{r_{E}}=-\frac{GM}{2a}}$.

Using $\boldsymbol{v_{esc}=\sqrt{2GM/r_{E}}}$ and $\boldsymbol{u_{0}=v_{0}/v_{esc}}$, we obtain

$\large \boldsymbol{\frac{2a}{r_{E}}=\frac{1}{1-u_{0}^{2}}}$.

(3)

Figure 5a shows the missile trajectory as an ellipse with focus $\boldsymbol{F_{1}}$ at the Earth’s center (the center of force).  The launch point P and the target point Q both lie on the planet’s surface, and the actual missile trajectory is the segment of the ellipse lying outside of the planet ($\boldsymbol{r\geq r_{E}}$).  The range of the missile is $\boldsymbol{R=2r_{E}\theta }$, where θ is the angle between $\boldsymbol{\overline{F_{1}P} }$ (or $\boldsymbol{\overline{F_{1}Q} }$) and the major axis of the ellipse.  The location of the second focus $\boldsymbol{F_{2}}$ depends on $\boldsymbol{\varepsilon }$, which is governed by the missile launch angle.  For all points P (and, by symmetry,  Q), $\boldsymbol{\overline{F_{1}P}+\overline{F_{2}P}=2a}$, regardless of the location of $\boldsymbol{F_{2}}$ along the major axis.  Since $\boldsymbol{\overline{F_{1}P}=r_{E}}$, then $\boldsymbol{\overline{F_{2}P}=2a-r_{E}}$.  By inspection, the angle θ is maximized when $\boldsymbol{\overline{F_{2}P}}$ is perpendicular to the major axis.  Therefore,

$\large \boldsymbol{sin\theta _{max}=\frac{2a-r_{E}}{r_{E}}=\frac{2a}{r_{E}}-1=\frac{u_{0}^{2}}{1-u{_{0}^{2}}}}$

and

$\large \boldsymbol{R_{max}=2r_{E}\theta _{max}}$ .

(4)

Besides calculating the maximum range, it is interesting to find the launch angle associated with it.  In Figure 5b, let $\boldsymbol{\overline{mn}}$ once again be tangent to the ellipse at the launch point P.  Using the reflection property discussed above, angle β = β′, so $\boldsymbol{\alpha +2\beta =\pi }$.  (Note: β is the launch angle measured with respect to the local vertical direction at P.)   Since  $\boldsymbol{F_{1}PF_{2}}$ is a right triangle, angle $\boldsymbol{\alpha =\pi /2 -\theta _{max}}$.  Therefore,

$\large \boldsymbol{\beta =\frac{\pi }{2}-\left ( \frac{\pi }{4}-\frac{\theta _{max}}{2} \right )=\frac{\pi }{4}+\frac{\theta _{max}}{2}}$ .

In the literature, the launch angle is more commonly measured relative to the local horizontal, i.e., $\boldsymbol{\varphi =\pi /2 - \beta }$, or

$\large \boldsymbol{\varphi =\frac{\pi }{4}-\frac{\theta _{max}}{2}}$ .

(5)

Note that as $\boldsymbol{u_{0}\rightarrow 0}$$\boldsymbol{\theta _{max}\rightarrow 0}$ and $\boldsymbol{\varphi \rightarrow \pi /4}$, or 45°, just as we would expect from basic kinematics.  If $\boldsymbol{u_{0}=1/\sqrt{2}}$ (the condition for low Earth orbit), $\boldsymbol{\theta _{max}=\pi /2}$ and $\boldsymbol{\varphi =0}$, i.e., the missile should be launched parallel to the Earth’s surface.

According to US Pacific Command, the Hwasong-14 missile launched by North Korea on July 4 was aloft for 37 minutes.  From Figure 2, this indicates that $\boldsymbol{u_{0}=0.57}$, or $\boldsymbol{v_{0}=6.4}$ km/s.  Eqn. 1 then tells us that the missile reached a maximum altitude $\boldsymbol{h_{max}\simeq 3000}$ km during its test flight.  From Eqn. 4, $\boldsymbol{sin\theta _{max}=0.48}$, or $\boldsymbol{\theta _{max}=0.50}$ (radians), so the maximum range is $\boldsymbol{R_{max}=2r_{E}\theta _{max}\simeq 6400}$ km.  The launch angle associated with this range is $\boldsymbol{\varphi =\pi /4 - 0.25 = 0.54}$ rad, or 31°.  These results are in good agreement with those of Dr. David Wright, Senior Scientist at the Union of Concerned Scientists, who estimated (Ref. 7) the maximum altitude to be “more than 2800 km” and the maximum range to be “roughly 6700 km.”

### D.  Questions for Students

1. The July 4 missile flew for 37 minutes and splashed down about 950 km from its launch site.  Show that the launch angle must have been close to 90°.  (Ans: $\boldsymbol{\bar{v_{\theta }}=950/(37\times 60)=0.43}$ km/s, while $\boldsymbol{v_{0}=6.4}$ km/s.  $\boldsymbol{cos\varphi =0.43/6.4 = 0.067}$, or $\boldsymbol{\varphi \simeq 86}$°)

2.   The July 28 missile flew for 47 minutes.  Find the maximum altitude during the test flight, the maximum range  $\boldsymbol{R_{max}}$, and the optimum launch angle $\boldsymbol{\varphi }$.  (Ans: $\boldsymbol{h_{max}\simeq 4200}$ km, $\boldsymbol{R_{max}\simeq 9200}$ km, and $\boldsymbol{\varphi \simeq 24}$°)

3.  Your calculation of $\boldsymbol{R_{max}}$ above assumed that burnout occurred at zero altitude.  Realistically, burnout occurs at an altitude of several hundred kilometers.  The “official” estimate of the range of the July 28 missile was 10,400 km.  Using the launch angle calculated above, estimate the burnout altitude of the missile.  (Ans: about 250 km)

### References

1. Choe Sang-Hun, “U.S. Confirms North Korea Fires an Intercontinental Ballistic Missile,” New York Times, July 4, 2017
2. Choe Sang-Hun and Eileen Sullivan, “North Korea Launches Ballistic Missile, the Pentagon Says,” New York Times, July 28, 2017
3. Robert W. Stanley, “Numerical methods in mechanics,” Am. J. Phys. 52(6), 499-507 (1984)
4. S. K. Foong, “From Moon-falls to motions under inverse square laws,” Eur. J. Phys. 29, 987-1003 (2008)
5. For a simple proof of the reflection property, see Physics from Planet Earth, J. C. Amato and E. J. Galvez, Taylor & Francis/CRC Press (2015) p. 471
6. ibid, section 8.2, “Kepler’s Laws”
7.  David Wright, Union of Concerned Scientists blog, July 3, 2017:  http://allthingsnuclear.org/dwright/north-korea-appears-to-launch-missile-with-6700-km-range

### Appendix

The series approximation shown in Figures 2 and 3 is derived by Foong in Ref. 4.  If the missile is at its peak altitude at time $\boldsymbol{t=0}$, then $\boldsymbol{r(t)=r(-t)}$; i.e., $\boldsymbol{r(t)}$ is an even function of time, and may be expressed as a power series:

$\large \boldsymbol{r(t)=r_{max}+c_{2}t^{2}+c_{4}t^{4}+\cdots }$ .

The coefficients $\boldsymbol{c_{2n} }$ are chosen to satisfy the equation of motion

$\large \boldsymbol{\frac{d^{2}r}{dt^{2}}=-\frac{GM}{r^{2}}}$ .

Following Foong, we introduce dimensionless variables $\boldsymbol{x=r/r_{max}}$ and $\boldsymbol{\tau =t/T}$, where $\boldsymbol{T=\sqrt{r_{max}^{3}/GM}}$.  We then rewrite the power series as $\boldsymbol{x=1+{c_{2}}'\tau ^{2}+{c_{4}}'\tau ^{4}+\cdots }$, which must satisfy the equation of motion expressed in dimensionless units:

$\large \boldsymbol{1+x^{2}\frac{d^{2}x}{d\tau ^{2}}=0}$

or

$\large \boldsymbol{1+(1+{c_{2}}'\tau ^{2}+{c_{4}}'\tau ^{4}+\cdots )^{2}(2{c_{2}}'+12{c_{4}}'\tau ^{2}+\cdots )=0}$ .

Expanding the parentheses and multiplying,

$\large \boldsymbol{\left ( 1+2{c_{2}}' \right )+ \left ( 4{c_{2}}'^{2} + 12{c_{4}}' \right )\tau ^{2} +\cdots =0}$ .

For the equality to hold over the entire flight time of the missile, the coefficient of each power of τ must equal zero, so $\boldsymbol{{c_{2}}'=-1/2}$ and $\boldsymbol{{c_{4}}'=-1/12}$.  See Ref. 4 for the values of higher order coefficients.  Returning to the physical variables r and t, our power series solution is

$\large \boldsymbol{r=r_{max}-\frac{1}{2}\left ( \frac{GM}{r_{max}^{2}} \right )t^{2}-\frac{1}{12}\frac{1}{r_{max}^{3}}\left ( \frac{GM}{r_{max}^{2}} \right )^{2}t^{4}+\cdots }$

where $\boldsymbol{r_{max}}$ is given as a function of $\boldsymbol{u_{0}}$ in Eqn. 3.  The truncated series containing just the terms shown above is plotted in Fig. 2 for $\boldsymbol{u_{0}}$ = 0.5, 0.6 and 0.7.  The agreement with the numerical solution is quite good for $\boldsymbol{u_{0}< 0.5}$.

# The Trappist-1 Solar System: Seven Earth-sized Planets Potentially Harboring Life

## A. Introduction

The recent discovery (Ref. 1) of seven exoplanets orbiting the red dwarf star Trappist-1 is an exciting and potentially key development in the search for extraterrestrial life in the universe.  All of the seven planets are Earth-sized, and at least three lie in the temperate “Goldilocks” zone, neither too near nor too far from their host star for surface oceans to exist.  Given the mass and temperature of these planets, their atmospheres are likely to contain nitrogen, oxygen, water vapor, and carbon dioxide, just like on Earth.  As an added bonus, the small size of the host star will help astronomers determine the chemical composition of each planet’s atmosphere.

The goal of this post is to show students how scientists use elementary physics and trigonometry to derive many of the physical properties of this and other extrasolar systems.  We present straightforward explanations appropriate to an introductory physics course, and offer exercises suitable for a homework assignment or recitation class to help students appreciate the significance of the Trappist-1 discovery.  We hope these materials will be useful to you, and that your students will enjoy and benefit from them.

## B.  Background

In February 2017, an international team of astronomers, led by Michaël Gillon of Belgium, announced the discovery (Ref. 1) of seven Earth-sized exoplanets orbiting the red dwarf Trappist-1, so named for the instrument in Chile (TRAPPIST: Transiting Planets and Planetesimals Small Telescope) first employed in their work.  Observations continued for one year using an array of ground-based telescopes located in Chile, Morocco, Hawaii, Spain, and South Africa, plus the orbiting Spitzer Space Telescope.  The exoplanet discoveries were made by the transit photometry technique:  when a planet eclipses its host star, i.e., as it passes between the star and Earth, it blocks a tiny fraction of starlight from reaching Earth.  The resulting light curve (flux vs. time, see Fig. 2) and its period P (time between successive transits), plus estimates of the star’s mass $\boldsymbol{M_{\ast }}$ and radius $\boldsymbol{R_{\ast }}$, suffice to determine the planet’s radius $\boldsymbol{R_{P}}$ and the radius of its orbit  $\boldsymbol{a}$.  This information can then be used to estimate the surface temperature $\boldsymbol{T}$ and the potential for life on the planet.

Astronomers determine stellar masses by studying binary stars.  For each binary system, they measure the distance from Earth by parallax, and measure the apparent brightness of each star to calculate its luminosity $\boldsymbol{L_{\ast }}$.  Mass is found by studying the orbital motion of the binary, as described in PPE, Section 8.11.  By compiling data on many binaries, an empirical mass-luminosity relation $\boldsymbol{L_{\ast }=f(M_{\ast })}$ can be constructed.  A mass-radius relation can be constructed in a similar way, using data from eclipsing binary stars and long baseline interferometry.  For now, let’s accept the mass and radius of Trappist-1 given in Ref. 1: $\boldsymbol{M_{\ast }/M_{\bigodot }=0.0802\pm .0073}$$\boldsymbol{R_{\ast }/R_{\bigodot }=0.117\pm 0.0036}$, where the subscript $\boldsymbol{\bigodot }$ indicates the Sun.  Table 1 is adopted from Ref. 1, and lists many of the properties of the Trappist-1 exosystem.

1.   A commonly quoted mass-luminosity relation, suitable for some small-mass main sequence stars, is $\boldsymbol{L_{\ast }/L_{\bigodot }=(M_{\ast }/M_{\bigodot })^{2.62}}$, where the subscript $\boldsymbol{\bigodot }$ indicates the Sun.  The measured luminosity of Trappist-1 is $\boldsymbol{L_{\ast }/L_{\bigodot }=5.24\times 10^{-4}}$.  Show that this leads to a value of $\boldsymbol{M_{\ast }}$ that agrees with the value quoted in Ref. 1.

## C.  Kepler’s Laws and Orbital Radii

Assume that the mass of each planet is much less than that of the host star, $\boldsymbol{M_{p}\ll M_{\ast }}$, and that each planet is in a near-circular orbit of radius $\boldsymbol{a\gg R_{\ast }}$, just like the planets of our own solar system.  Kepler’s 3rd law is

$\LARGE \boldsymbol{\frac{P^{2}}{a^{3}}=\frac{4\pi ^{2}}{GM_{\ast }}}$

(1)

where P is the orbital period.  After calculating $\boldsymbol{M_{\ast }}$ from its luminosity, and measuring P , we can use Eqn. 1 to derive a.

2.  Astronomers have recently (February 2017) discovered seven exoplanets orbiting the star Trappist-1, an ultracool red dwarf about 12 pc from Earth.  This is an exciting discovery because at least three  planets are in the “temperate zone” of the host star, with surface temperatures favorable for life.  The table below lists properties of six bodies that orbit either Trappist-1 ($\boldsymbol{M_{\ast }\simeq 0.08M_{\bigodot }}$) or one other central body with mass $\boldsymbol{M\simeq 0.001M_{\bigodot }}$.  Fill in the missing numbers.  (Hint: Use ratios rather than converting to meters, seconds and kilograms.)  Can you identify M and the bodies in orbit about it?

## D.  Analyzing the Light Curve

Fig. 2 contains a great deal of information that can be used to further characterize the Trappist-1 system.  Even if $\boldsymbol{M_{\ast }}$ and $\boldsymbol{R_{\ast }}$ are not known a priori, the period P and the shape of the transiting planet’s light curve provide enough information to determine $\boldsymbol{R_{p}}$, the orbital radius a and inclination i, and the density of the host star $\boldsymbol{\rho _{\ast }}$.  Combining the star’s density with the appropriate mass-radius relation for red dwarfs yields $\boldsymbol{M_{\ast }}$ and $\boldsymbol{R_{\ast }}$.  This strategy uses Kepler’s 3rd law plus simple kinematics and trigonometry, making it amenable to first-year physics students.  It is described in detail by Seager and Mallén-Ornelas in Ref. 2, and we will follow their treatment closely.

As usual, assume that the exoplanet’s mass is much less than the star’s, and that it is in a circular orbit of radius $\boldsymbol{a\gg R_{\ast }}$.  Fig. 3a is a “side view” of the orbit, perpendicular to the line of sight from Earth to the star.  The inclination angle is defined so that if $\boldsymbol{i=90\: \textup{deg}}$, observers on Earth view the orbit edge-on.  The impact parameter b is defined by $\boldsymbol{bR_{\ast }=acosi}$.  Clearly, $\boldsymbol{b<1}$ for a transit to occur, and since $\boldsymbol{a\gg R_{\ast }}$, this implies that $\boldsymbol{i\simeq 90\; \textup{deg}}$.

From Fig. 3b, the amount of light blocked by the planet is proportional to its area, so the fractional change in flux reaching Earth during a transit is

$\LARGE \boldsymbol{\frac{\Delta F}{F}=\frac{R_{p}^{2}}{R_{\ast }^{2}}}$.

(2)

In practice, the straight line segments of the light curve shown in Fig. 3b are rounded by limb-darkening (as seen in Fig. 2): less light reaches us from areas near the edge of the stellar disc than from equal areas near the center of the disc.  For simplicity, we will ignore this effect.

3.  Fig. 2 displays the light curves of the seven planets orbiting Trappist-1.  Which planet has the largest radius?  (Ans: Trappist-1b)

Since $\boldsymbol{a\gg R_{\ast }}$, we may approximate the trajectory of the planet as a straight line as it moves across the face of the star.  For a circular orbit, the planet’s speed is $\boldsymbol{v=2\pi a/P}$.  By the Pythagorean theorem, the distance traveled by the planet in time $\boldsymbol{t_{t}}$ is $\boldsymbol{d_{t}=vt_{t}=2\sqrt{(R_{\ast }+R_{p})^{2}-b^{2}R_{\ast }^{2}}}$, so

$\LARGE \boldsymbol{t_{t}=\frac{P}{\pi a}\sqrt{(R_{\ast }+R_{p})^{2}-b^{2}R_{\ast }^{2}}}$.

(3)

4.  Derive the corresponding equation for $\large \boldsymbol{t_{f}}$
(Ans: $\large \boldsymbol{t_{t}=\frac{P}{\pi a}\sqrt{(R_{\ast }-R_{p})^{2}-b^{2}R_{\ast }^{2}}}$)

If $\boldsymbol{R_{\ast }}$ and $\boldsymbol{M_{\ast }}$ are not known a priori, then the only information available is the light curve and its period.  Nevertheless, it is still possible to proceed.  From Eqn. 3 and the corresponding expression for $\boldsymbol{t_{f}}$ (Question 4), it is easy to show that

$\LARGE \boldsymbol{(t_{t}^{2}-t_{f}^{2})^{\frac{1}{2}}=\frac{2P}{\pi a}R_{\ast }\left ( \frac{\Delta F}{F} \right )^{\frac{1}{4}}}$,

where we have used Eqn. 2 to express $\boldsymbol{R_{p}/R_{\ast }=\sqrt{\Delta F/F}}$.  Rearranging the above expression,

$\LARGE\boldsymbol{\frac{a}{R_{\ast }}=\frac{2P}{\pi }\frac{(\Delta F/F)^{\frac{1}{4}}}{(t_{t}^{2}-t_{f}^{2})^{\frac{1}{2}}}}$

or

$\LARGE \boldsymbol{\frac{a^{3}}{R_{\ast }^{3}}=\frac{8P^{3}}{\pi^{3} }\frac{(\Delta F/F)^{\frac{3}{4}}}{(t_{t}^{2}-t_{f}^{2})^{\frac{3}{2}}}}$.

(4)

Combining Eqn. 4 with $\boldsymbol{M_{\ast }=4\pi ^{2}a^{3}/GP^{2}}$ (Kepler’s 3rd law) yields an expression for the stellar “density”:

$\LARGE \boldsymbol{\rho _{\ast }=\frac{M_{\ast }}{R_{\ast }^{3}}=\frac{32P}{\pi G}\frac{(\Delta F/F)^{\frac{3}{4}}}{(t_{t}^{2}-t_{f}^{2})^{\frac{3}{2}}}}$.

It is convenient to re-express this in terms of the solar “density”

$\LARGE \boldsymbol{\frac{\rho _{\ast }}{\rho _{\bigodot }}=\frac{32R_{\bigodot }^{3}}{\pi GM_{\bigodot }}\: \frac{(\Delta F/F)^{\frac{3}{4}}}{(t_{t}^{2}-t_{f}^{2})^{\frac{3}{2}}}\: P}$

(5)

where the pre-factor $\boldsymbol{32R_{\bigodot }^{3}/\pi GM_{\bigodot }=3.46\times 10^{-3}\: {\textrm{day}}^{2}}$.  This is the key result of this section.

For classroom purposes, we can crudely express the mass-radius relation as $\boldsymbol{R_{\ast }=M_{\ast }^{x}}$. Choosing $\boldsymbol{x=0.85}$, we can now calculate the mass and radius of Trappist-1 in solar mass units:

$\LARGE\boldsymbol{\frac{M_{\ast }}{M_{\bigodot }}=\frac{\rho _{\ast }}{\rho _{\bigodot }}\frac{R_{\ast }^{3}}{R_{\bigodot }^{3}}=\frac{\rho _{\ast }}{\rho _{\bigodot }}\left ( \frac{M_{\ast }}{M_{\bigodot }} \right )^{2.55}}$

or

$\LARGE \boldsymbol{\frac{M_{\ast }}{M_{\bigodot }}=\left ( \frac{\rho _{\ast }}{\rho _{\bigodot }} \right )^{-0.65}\; and\; \; \; \frac{R_{\ast }}{R_{\bigodot }}=\left ( \frac{\rho_{\ast } }{\rho_{\bigodot } } \right )^{-0.55}}$,

(6)

where $\boldsymbol{\rho _{\bigodot }\equiv M_{\bigodot }/R_{\bigodot }^{3}=5.90\times 10^{3}\: \mathrm{kg/m}^{3}}$.  The validity of these results depends sensitively on the mass-radius relation chosen.  Admittedly, the exponent $\boldsymbol{x=0.85}$ was chosen to yield results consistent with those of Ref. 1.  In the literature, $\boldsymbol{x=0.80}$ is more commonly quoted.

5.  According to Ref. 1, the transit duration of Trappist-1d is $\boldsymbol{t_{t}=49.1\: \textrm{min}}$, the period is $\boldsymbol{P=4.05\: \textrm{days}}$, and the transit depth is $\boldsymbol{\Delta F/F=0.0037}$.

a.  Use the values of $\boldsymbol{M_{\ast }}$ and $\boldsymbol{R_{\ast }}$ given in section B to compute $\boldsymbol{\rho _{\ast }/\rho _{\bigodot }}$ and then find $\boldsymbol{t_{f}}$.  Check your answer against the light curve in Fig. 2.  (Ans: $\boldsymbol{t_{f}\approx 0.88t_{t}=43.3\: \textrm{min}}$)

b.  Find the radius of the planet’s orbit.  (Ans: $\boldsymbol{a=2.14\times 10^{-2}\: \textup{\textrm{AU}}}$)

c.  Find the radius of the planet.  Express $\boldsymbol{R_{p}}$ as a fraction of $\boldsymbol{R_{Earth}}$.  (Ans: $\boldsymbol{R_{p}=0.78\, R_{Earth}}$)

6.  Researchers at Trump University recently reported the discovery of an eighth planet orbiting Trappist-1.  In a series of tweets, team leaders provided a summary of their observations taken with the privately funded Andrew Jackson Telescope in Palm Beach, FL:  $\boldsymbol{P=15.0\pm 0.2\: \textrm{min}}$$\boldsymbol{\Delta F/F=0.005\pm 0.0001}$$\boldsymbol{t_{t}=86.4\pm 0.2\: \textrm{min}}$, and $\boldsymbol{t_{f}=55.4\pm 0.2\: \textrm{min}}$.  Their measurements have not yet been confirmed.  If you were a program director at the National Science Foundation, would you approve a grant to the Trump exoplanet program?

## E.  Too Hot, Too Cold, or Just Right?

A blackbody is an object that absorbs all of the radiation incident on its surface.  (Since it reflects no light, it appears black.)  But a blackbody is not just an ideal absorber of radiant energy, it is an ideal emitter of energy as well.  The power S emitted by a blackbody per unit surface area is given by the Stefan-Boltzmann Law,

$\LARGE \boldsymbol{S=\sigma _{B}T^{4}}$,

(7)

where T is the surface temperature in kelvin, and $\boldsymbol{\sigma _{B}=5.67\times 10^{-8}\: \textrm{W/m}^{2}\textrm{-K}^{4}}$ is called the Stefan-Boltzmann constant.  The Sun and a cloudless Earth behave approximately as blackbodies.  The power emitted by the Sun, its luminosity, is $\boldsymbol{L_{\bigodot }=4\pi \sigma _{B}R_{\bigodot }^{2}T_{\bigodot }^{4}=3.83\times 10^{26}\: \textrm{W}}$.  At a distance r from the Sun, this power is spread evenly over a sphere of area $\boldsymbol{4\pi r^{2}}$.  A spherical planet of radius $\boldsymbol{R_{P}}$ presents a cross-sectional area $\boldsymbol{\pi R_{P}^{2}}$ to the Sun, so if it acts as a blackbody, it intercepts and absorbs energy at a rate $\boldsymbol{W_{abs}=L_{\bigodot }R_{P}^{2}/4r^{2}}$.  This energy warms the planet’s surface to a temperature , causing it to emit energy at a rate $\boldsymbol{W_{emit}=4\pi R_{P}^{2}\cdot \sigma _{B}T^{4}}$. When $\boldsymbol{W_{emit}=W_{abs}}$, the temperature reaches a steady-state value given by

$\LARGE \boldsymbol{T=\left ( \frac{L_{\bigodot }}{16\pi \sigma _{B}r^{2}} \right )^{\frac{1}{4}}}$.

(8)

7.  The Earth is 1 AU ($\boldsymbol{1.5\times 10^{11}\: \textrm{m}}$) from the Sun.  What is the steady-state temperature on Earth, according to Eqn. 8?  (Ans: 278 K.  Keep in mind, though, that this ignores the greenhouse effect, which raises the temperature significantly.)

To apply Eqn. 8 to an extrasolar system, simply replace $\boldsymbol{L_{\bigodot }}$ by the luminosity of the host star $\boldsymbol{L_{\ast }}$.

8.  The period of Trappist-1d is 4.05 days.  Use the star’s mass $\boldsymbol{M_{\ast }/M_{\bigodot }=0.08}$ and luminosity $\boldsymbol{L_{\ast }/L_{\bigodot }=5.24\times 10^{-4}}$ to calculate the steady-state surface temperature on the planet.  Based only on your answer, what is the likelihood of finding liquid water on this planet?  (Ans: T ≈ 300 K.  However, if the planet has an Earth-like atmosphere, the greenhouse effect would raise its temperature too high to harbor liquid water on its surface.  This is mentioned briefly in Ref. 1.)

9.  What is the likelihood of finding intelligent life on Trappist-1i, the planet claimed by researchers at Trump University?

## References

1.   Michaël Gillon et al, Nature 542, 456 (2017)

2.  S. Seager and G. Mallén-Ornelas, Ap. J. 585, 1038 (2003)

# Work, Energy and the “Satellite Drag Paradox”

## A.  Introduction

Here is a lovely application of energy conservation and the “work-energy theorem” that will surely intrigue your first-year physics students.  A satellite in low Earth orbit experiences a drag force $\boldsymbol{\vec{F}_{d}}$ as it passes through the upper  reaches of the planet’s atmosphere.  (See PPE, Section 6.10, p. 188 ff.)  The drag force is antiparallel to the satellite’s velocity $\boldsymbol{\vec{v}}$, so it robs the spacecraft of energy and causes it to spiral inward toward the Earth’s surface.  However, assuming a gradual decay and a near-circular trajectory, it is easy to show (see below) that the speed of the spacecraft must increase as the orbital radius decays.  It seems unphysical that by applying a drag (friction) force to a body, its speed will increase!

This “paradox” has been around since the 1950’s, but to my knowledge it is not treated in introductory physics textbooks.  One of the clearest and most complete solutions was given by B. D. Mills in 1959, and we will use his approach in the following paragraphs.  As usual, the presentation will be simplified as much as possible to make it tractable for introductory students.

## B.  The Satellite Drag Paradox

Let the satellite be moving on a low Earth near-circular trajectory with a slowly decreasing orbital radius $\boldsymbol{r(t)}$ due to the drag force $\boldsymbol{\vec{F}_{d}}$.  Mills (Ref. 1) shows that for a weak drag force, the tangential speed of the spacecraft is at all times nearly equal to that of a body moving in a perfectly circular orbit of the same instantaneous radius.  (See also Ref. 5.)  Let $\boldsymbol{M}$ and $\boldsymbol{m}$ denote the mass of the Earth and satellite, respectively.  From Newton’s Second Law, $\boldsymbol{mv^{2}/r=GMm/r}$, the tangential speed $\boldsymbol{v_{c}}$ is

$\large \boldsymbol{v_{c}=\sqrt{GM/r}}$

(1)

and it increases as the orbital radius decreases.  The kinetic energy of a body $\large \boldsymbol{m}$ in circular orbit is found using Eqn. 1:

$\large\boldsymbol{K=\frac{1}{2}mv_{c}^{2}=\frac{1}{2}m\frac{GM}{r}=-\frac{1}{2}U}$,

(2)

where $\boldsymbol{U=-GMm/r}$ is the potential energy, so the total energy of the orbiting body is

$\large \boldsymbol{E=K+U=\frac{1}{2}U=-\frac{GMm}{2r}}$.

(3)

These results should be familiar to students.

1.  Using Eqns. 2 and 3, verify that as $\boldsymbol{r}$ decreases, $\boldsymbol{K}$ increases even though $\boldsymbol{U}$ and $\boldsymbol{E}$ decrease.

If the spacecraft’s orbit is initially elliptical (as opposed to circular), it will soon circularize due to the atmospheric drag force.  The drag force increases dramatically at lower altitudes, so it will be greatest when the spacecraft passes through perigee. Imagine that the drag force produces an instantaneous negative impulse at perigee, with change in velocity $\boldsymbol{\Delta v}$ antiparallel to the spacecraft’s velocity.  This reduces the total energy of the spacecraft, shrinking its semi-major axis without changing its perigee position $\boldsymbol{r_{p}}$.  As a result, the apogee distance $\boldsymbol{r_{a}}$ shrinks, and as $\boldsymbol{r_{a}\rightarrow r_{p}}$, the orbit becomes more circular.  So even if the initial orbit is significantly elliptical, atmospheric drag will soon circularize it, and the analysis presented below will be applicable.

For a spacecraft experiencing a drag force $\boldsymbol{\vec{F}_{d}}$, the rate of change of $\boldsymbol{E}$ is equal to $\boldsymbol{\vec{F}_{d}\cdot \vec{v}}$, where $\boldsymbol{\vec{v}}$ is the total velocity of the spacecraft, including the inspiral motion.  This follows because $\boldsymbol{\vec{F}_{d}}$ is the only external force acting on the Earth-satellite system.  Taking the time derivative of Eqn. 3, and noting that $\boldsymbol{\vec{F}_{d}}$ is antiparallel to the velocity $\boldsymbol{\vec{v}}$, we find

$\large \boldsymbol{\frac{dE}{dt}=\frac{GMm}{2r^{2}}\frac{dr}{dt}=-F_{d}v}$

or

$\large \boldsymbol{\frac{dr}{dt}=-\frac{2r^{2}}{GMm}F_{d}v}$.

Figure 2 illustrates the forces acting on the satellite, and the components of its velocity: $\boldsymbol{\vec{v}=v_{c}\hat{i}+\frac{dr}{dt}\hat{j}}$.  For satellites passing through the thin upper atmosphere of a planet, $\boldsymbol{\left | \frac{dr}{dt} \right |\ll v_{c}}$, so we can approximate $\boldsymbol{\vec{F}_{d}}$ to be antiparallel to $\boldsymbol{\vec{v}_{c}}$, as shown.  Thus, the above result may be written as

$\large \boldsymbol{\frac{dr}{dt}=-\frac{2r^{2}}{GMm}F_{d}v_{c}}$.

(4)

Using the work-energy theorem $\boldsymbol{\Delta K=\vec{F}\cdot \Delta \vec{r}}$, where $\boldsymbol{\vec{F}}$ is the total force acting on the body, we obtain

$\large \boldsymbol{\frac{dK}{dt}=\vec{F}\cdot \vec{v}=F_{g}\frac{dr}{dt}-F_{d}v_{c}}$

(5)

where $\boldsymbol{F_{g}=-GMm/r^{2}}$ is the gravitational force exerted on the spacecraft by the Earth.  Eqn. 5 reveals the source of the apparent paradox: while the drag force decreases $\boldsymbol{K}$, the satellite gains speed as it falls toward Earth.  To evaluate $\boldsymbol{dK/dt}$, substitute Eqn. 4 into Eqn. 5:

$\large \boldsymbol{\frac{dK}{dt}=-\frac{GMm}{r^{2}}\left ( -\frac{2r^{2}}{GMm}F_{d}v_{c} \right )-F_{d}v_{c}=F_{d}v_{c}}$.

(6)

Surprisingly, the drag term appears in Eqn. 6 without a minus sign!  As described by Mills, it is “as if the air drag force, reversed, were pushing the satellite.”

## C.  Application: The International Space Station (ISS)

The International Space Station (Fig. 1) is maintained in a near-circular orbit of approximate altitude 400 km.  Figure 3 (Ref. 2) is a plot of the ISS’s mean altitude vs. time, and shows how the altitude is periodically boosted (the near-vertical line segments) and how it decays between boosts.  For an on-board demonstration of what ISS astronauts experience during a “reboost,” see Ref. 4.

In November 2016, the ISS was boosted to a mean altitude of 406.5 km, and then fell steadily to 404.5 km over the subsequent three month interval (92 d).

3.  What is the orbital speed of the ISS at an altitude of about 400 km?  (Ans: 7.66 km/s)

4.  From November 2016 to February 2017, the altitude decreased by 2 km.  Verify that $\boldsymbol{\left | dr/dt \right |\ll v_{c}}$, as asserted above.

5.  Calculate the change in the orbital speed $\boldsymbol{v_{c}}$ during that three month period.  (Hint:  $\boldsymbol{\Delta v/v_{c}=(?)\Delta r/r}$.  Ans: $\boldsymbol{\Delta v_{c}=1.13\: \textup{m/s}}$)

Now let’s calculate how much rocket fuel must be consumed to boost the ISS back to its original altitude of 406.5 km.  To keep things simple, let’s assume that the boost is done in a single short burn, as shown in Figure 4.  A single burn will change a circular orbit to an elliptical one, with perigee $\boldsymbol{r_{p}=R_{E}+404.5 \: \textup{km}}$ at the site of the burn.  If we wish the semi-major axis to be $\boldsymbol{a=R_{E}+406.5\: \textup{km}}$, then the apogee must be $\boldsymbol{r_{a}=R_{E}+408.5\: \textup{km}}$, located 180º from the site of the burn.  This results in an eccentricity $\boldsymbol{\epsilon \approx 3\times 10^{-4}}$, which we will ignore.  (To improve the circularity of the boosted orbit, a two burn strategy is sometimes used.  See Ref. 3.)

The change in the ISS’s orbital energy due to the burn is

$\large \boldsymbol{E_{f}-E_{i}=\frac{1}{2}m(v_{i}+\Delta v)^{2}-\frac{1}{2}mv_{i}^{2}\simeq mv_{i}\Delta v}$,

where $\boldsymbol{v_{i}}$ is the speed before the burn, and $\boldsymbol{v_{i}+\Delta v}$ is the speed afterwards.  We can solve for $\boldsymbol{\Delta v}$ using Eqn. 3, or, since the altitude change is so small, set $\boldsymbol{\Delta U=mg\Delta a}$, where $\boldsymbol{g}$ is the gravitational acceleration at the altitude of the ISS, and $\boldsymbol{\Delta a=406.5-404.5=2.0\: \textup{km}}$ is the change in the mean altitude.  From Eqn. 3, $\boldsymbol{\Delta E=\frac{1}{2}\Delta U}$.

6.  What is the change of velocity needed to boost the altitude of the ISS by 2.0 km?  (Ans: 1.13 km/s)

ISS orbit boosts are carried out using thrusters on the cargo spacecraft that supply the ISS, such as the Russian Progress modules.  (In the past, the US Space Shuttles were used.)  The Progress thrusters use a fuel (UDMH + NTO) with exhaust velocity 2.7 km/s.

7.  The mass of the ISS is $\boldsymbol{4.2\times 10^{5}\: \textup{kg}}$.  Calculate the fuel mass needed to boost its altitude by 2.0 km.  (Ans: 176 kg)

8.  Calculate the magnitude of the drag force $F_{d}$ acting on the ISS during November 2016 to February 2017.  You might be surprised by your answer!  (Ans: .06 N)

9.  In ref. 4, astronaut Jeff Williams discusses what happens during a typical orbit boost.  For the maneuver shown in the video, $\boldsymbol{\Delta v=2.7\: \textup{m/s}}$, and the acceleration was $\boldsymbol{1.85\: cm/s^{2}}$.  What was the duration of the burn?  (Ans: 146 s)

10.  Does your answer to Question 7 include the fuel required to overcome the drag that lowered the orbit in the first place?  Why is this a small correction?  (Ans: the duration of the reboost is only a few minutes, whereas the orbital period is about 90 minutes.  It took many orbits for the drag force to lower the ISS’s altitude.)

11.  There are at least two minor misstatements in the video of Ref. 4.  Can you spot them?  (Ans: At 1:09, Williams states that because of atmospheric drag, “over a period of time we slow down and our altitude over the Earth decreases.  At 3:04 he states “We’re in weightlessness right now, and there’s no acceleration …”  Actually, he and the ISS are undergoing centripetal acceleration of about $\boldsymbol{8.66\: m/s^{2}}$ towards the Earth.)

## Acknowledgment

After writing a draft of this post, I asked Dr. Philip Blanco – a valued contributor to this blog – to review it for correctness and to offer suggestions for improvement.  He not only did that, but also supplied a numerical simulation of the trajectory of an ISS-like spacecraft immersed in a thin atmosphere.  His solution, shown below (Figure 5), shows that the spacecraft indeed follows a near-circular orbit right up to the point where a catastrophe occurs, whence it plunges toward Earth.  It would be interesting to calculate when that catastrophe occurs.

My thanks to Dr. Blanco for his assistance.  Of course, I bear the full responsibility for any remaining errors.

## References

1.  Blake D. Mills, Am. J. Phys. 27, 115 (1959)

2.   Chris Peat, http://www.heavens-above.com/IssHeight.aspx

4.  https://www.youtube.com/watch?v=sI8ldDyr3G0, courtesy NASA Johnson Space Center

5.  Leon Blitzer, Am. J. Phys. 39, 882 (1971)