# 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)

# Rosetta’s Final Act: Free-Fall to the Surface of Comet 67P

The European Space Agency’s historic Rosetta mission is over.  On September 30, 2016, the spacecraft was allowed to fall ballistically (unpowered) to the surface of comet 67P/Churyumov-Gerasimenko, capturing numerous high-resolution photos of the comet’s surface during its 14 hour descent.  The final image, taken from an altitude of 20 meters, is shown in Figure 2.

Comets are presumed to be the primordial building blocks of the solar system, so by observing 67P at close range, Rosetta was investigating the origins of our solar system.

During its 12 year lifespan, Rosetta became the first spacecraft to orbit a comet and observe the chemical processes taking place as it passed through perihelion in August 2015, when the rate of gas and dust evolution peaked.  In November 2014, Rosetta was the first spacecraft to deploy a lander, Philae, to the surface of a comet.  Although Philae functioned only briefly, the two allied space probes made numerous discoveries about the chemical makeup and geology of comets.

Comet 67P has a period of 6.45 years and an aphelion distance of 5.7 AU.  At that distance, Rosetta‘s onboard heating system would have been incapable of keeping the spacecraft warm enough for its instruments to hibernate safely.  Consequently, flight directors opted to end the mission by crash-landing the spacecraft onto the comet, gathering as much data as possible while its instruments were still functioning.  Today, 67P is 4.3 AU from the Sun.  As it heads toward aphelion in November, 2018, it carries the remains of Rosetta and Philae, the two pioneering spacecraft that doggedly chased it for over a decade.

In this blog post, we present a set of student-ready exercises based on Rosetta‘s final maneuver and descent to the surface of 67P.  These exercises were inspired by Dr. Philip Blanco’s recent letter (Ref. 1) to the editor of The Physics Teacher.  Both Philip and I are deeply indebted to Pablo Muñoz, flight dynamics engineer at the ESA’s European Space Operations Centre in Darmstadt, Germany.  Pablo supplied indispensible data on the position and velocity of the spacecraft before and after its final thruster burn, plus the mass of Rosetta and the amount of fuel consumed during that final burn.  Without his input, none of these exercises could have been written.

To make the exercises suitable for introductory students, I have modified the details slightly, as indicated below.   To fully appreciate the complexity of Rosetta‘s dance with its cometary partner, follow the links in Ref. 4 and 5 to the beautiful ESA videos of the mission.  Have fun!

## A.  Problems

In these exercises, $\boldsymbol{ GM=666.15\mathrm{ m^{3}/s^{2}}}$, where
$\boldsymbol{ M}$ is the mass of the comet ($\boldsymbol{ M\simeq 1.0\times 10^{13}} \textup{ kg}$), and model the comet as a homogeneous sphere of radius 2.00 km.

Prior to its plunge to the comet surface, Rosetta was placed in an unpowered elliptical orbit with apocenter (farthest distance to the comet’s CM) $\boldsymbol{r_{a}=23.44}\textup{ km}$ and pericenter $\boldsymbol{r_{p}=13.58}\textup{ km}$.  This is shown nicely in Ref. 4.

1.  Calculate the specific energy $\boldsymbol{ E/m}$ (where $\boldsymbol{m}$ is Rosetta‘s mass) and the orbital period of the spacecraft.  (Ans: $\boldsymbol{E/m=-1.800\times 10^{-2}}\textup{ J/kg}$$\boldsymbol{ T=6.13\times 10^{5}}\textup{ s}\simeq 1\textup{ week}$)

2.  What was the speed $\boldsymbol{ v_{0}}$ of the spacecraft at apocenter? Compare this to the speed of a spacecraft in low Earth orbit.  Why the big difference?  (Ans: $\boldsymbol{ v_{0}=0.1444}\textup{ m/s}$$\boldsymbol{M_{67P}/M_{Earth}\sim 10^{-12}}$)

The final maneuver that sent Rosetta diving to the comet surface occurred at a true anomaly of 208º, i.e., 28º beyond apocenter (0º corresponds to pericenter).  (See Ref. 4.)  To streamline the following questions, assume instead that the maneuver took place exactly at apocenter.  This alters the results modestly.

3.   Immediately after the burn, Rosetta‘s velocity was $\boldsymbol{ \vec{v}{_{0}}'=-0.0211\hat{i}-0.3255\hat{j}}\textup{ m/s}$, where $\boldsymbol{ \hat{i}}$ and $\boldsymbol{ \hat{j}}$ refer to the radial and azimuthal directions (Ref. 2) at apocenter.  See Figure 3.  Show that the burn inserted the spacecraft into an unbounded hyperbolic trajectory.  (Ans: $\boldsymbol{ {E}'/m=2.478\times 10^{-2}}\textup{ J/kg}>0$)

4.  Note that $\boldsymbol{ {\vec{v}}'_{0}}$ is not perfectly radial, i.e., directed toward the comet’s CM.  Convince someone that the spacecraft will still strike the comet.  This question and the solution below were suggested by Philip Blanco.          (Ans:  Compare the angle subtended by the comet at apocenter to the angle between the velocity vector $\boldsymbol{{\vec{v}}'_{0}}$ and the radial direction.)

5.  With what speed $\boldsymbol{ v}$ will the spacecraft strike the surface?  Compare this to the surface escape speed.  (Ans: $\boldsymbol{ v=0.846}\textup{ m/s}$$\boldsymbol{ v_{esc}=0.816}\textup{ m/s}$)

6.  Find the spacecraft’s velocity $\boldsymbol{ \vec{v}=v_{x}\hat{i}+\vec{v}_{y}\hat{j}}$ just before it strikes the surface.  Hint: first conserve angular momentum.  (Ans: $\boldsymbol{ \vec{v}=-0.809\hat{i}-0.247\hat{j}}\textup{ m/s}$)

7.  Calculate the magnitude of the change in velocity caused by the thruster burn.  The mass of the spacecraft at this time was 1422.5 kg, and the burn consumed 0.2 kg of rocket fuel.  What was the effective exhaust velocity of the fuel (a mixture of monomethyl hydrazine (MMH) and nitrogen tetroxide ($\large \mathrm{N_{2}O_{4}}$ or NOX))?  (Ans: $\boldsymbol{ \Delta v_{0}=.365}\textup{ m/s}$$\boldsymbol{ v_{ex}=2.6}\textup{ km/s}$.  This is substantially less than the “published” value of 3.35 km/s (Ref. 3) for MMH/NOX because a fuel leak forced flight engineers to operate the thrusters at lower than optimal pressure.)

8.  In September 2016, comet 67P has a rotation period of 12.055 hours.  This is less than its period in 2014 (12.4 hours), when Philae was released, due to outgassing from the comet as it passed through perihelion.   Assume that the comet’s rotation axis was perpendicular to the plane of the spacecraft’s motion (Figure 3).  Rosetta‘s velocity $\boldsymbol{ \vec{v}_{0}{}'}$ was selected so that its impact velocity relative to the rotating surface was nearly vertical, i.e., perpendicular to the impact surface.  (Hence, an observer standing on the surface at the impact point would have seen the spacecraft descending vertically.)  Show that this is roughly consistent with your calculated value of $\boldsymbol{v_{x}}$ in Question 6.  (Note: because the comet has such an irregular shape, the local vertical direction was significantly different from the radial direction.) (Ans:$\boldsymbol{ v_{rot}=-\omega R=-29.0\textup{ cm/s} \approx v_{x}=-24.7\textup{ cm/s}}$)

The following question is suitable for users of Physics from Planet Earth who have studied Section 12.9: Application: the “Slingshot Effect” Revisited, or who have seen a comparable treatment elsewhere.

9.  Suppose Rosetta‘s final maneuver had been incorrectly executed, and the spacecraft just missed striking the comet, instead grazing its surface as it flew by.

a.  From Question 3, we know that Rosetta‘s speed after the burn was 0.3262 m/s at apocenter.  In the comet’s reference frame, what will be its speed after it has flown past the comet and is very far from it?  (Ans: $\boldsymbol{ v_{\infty }=0.2226\textup{ m/s}}$)

b.  By what angle Θ will the spacecraft be deflected after its close fly-by of the comet?  (Ans:  $\boldsymbol{ \tan (\Theta /2)=1.769}$$\boldsymbol{ \Theta =121^{\circ}}$.  See Eqn. 12.12 of PPE.)

## B.  Sphere of Influence of Comet 67P

Here is an additional “space-themed” topic that is quite suitable for discussion in an introductory mechanics course.  As far as I know, it is not discussed in any currently available textbook (including PPE).  I thank Philip Blanco for suggesting it.

So far, we have ignored the influence of the Sun on the motion of Rosetta about comet 67P.  We might justify this by invoking Einstein’s Principle of Equivalence (see Section 7.5 of PPE) which states that in a free-fall reference frame, the local effects of gravity disappear.  Comet 67P’s heliocentric orbit is free-fall motion about the Sun, so in the comet’s reference frame, we can ignore the Sun’s gravitational force on the comet (which is obvious).  However, there is a small but important complication when we consider the motion of the spacecraft in the same reference frame.  Since the Sun’s gravitational force varies inversely with the square of the distance, the spacecraft’s acceleration toward the Sun varies as the distance between it and the Sun changes along its orbit.  Let $\boldsymbol{ R}$ be the distance from the Sun to the comet, and $\boldsymbol{ r}$ be the distance from the comet’s CM to the spacecraft.  See Figure 4.  When the spacecraft is at position a(b) in the figure, its distance from the Sun is a maximum (minimum), and its acceleration due to the Sun is

$\LARGE \boldsymbol{a=\frac{GM_{Sun}}{(R\pm r)^{^{2}}}\simeq \frac{GM_{Sun}}{R^{2}}\left ( 1\mp \frac{2r}{R} \right )}$.

(1)

The acceleration of the comet toward the Sun is $\boldsymbol{ {a}'=GM_{Sun}/R^{2}}$, so the maximum difference between the spacecraft’s acceleration and that of the comet is $\boldsymbol{ \Delta a_{max}=2GM_{Sun}\, r/R^{3}}$.  How close is the spacecraft’s motion to an unperturbed Keplerian orbit about the comet?  The answer depends on the ratio between $\boldsymbol{ \Delta a_{max}}$ and the centripetal acceleration of the spacecraft toward the comet, $\boldsymbol{ a_{c}=GM/r^{2}}$, where $\boldsymbol{ M}$ is the comet mass.  The smaller this ratio is, the smaller the influence of the Sun, and the better the motion can be described as an unperturbed Keplerian orbit about the comet.

Now let’s ask the opposite question.  Imagine the spacecraft is far enough from the comet that its motion is better described as a heliocentric orbit perturbed slightly by the comet.  How close is its motion to a pure heliocentric orbit?  The answer depends on the ratio between the acceleration $\boldsymbol{ a_{c}}$ toward the comet and the acceleration toward the Sun $\boldsymbol{ {a}'}$:

$\LARGE\boldsymbol{ \frac{a_{c}}{{a}'}=\frac{M}{M_{Sun}}\frac{R^{2}}{r^{2}}}$.

(2)

The smaller this ratio is, the smaller the influence of the comet, and the better the motion can be described as a Keplerian orbit about the Sun.  The “sphere of influence” of the comet is the region within which the motion of the spacecraft is only weakly perturbed by the Sun.  By convention, its radius $\boldsymbol{ r_{soi}}$ is found by equating the ratios in Eqns. 1 and 2, and solving for $\boldsymbol{ r\textup{ }\: (=r_{soi})}$:

$\LARGE\boldsymbol{\frac{2M_{Sun}}{M}\frac{r^{3}}{R^{3}}=\frac{M}{M_{Sun}}\frac{R^{2}}{r^{2}}}$

or

$\LARGE\boldsymbol{ \frac{r^{5}}{R^{5}}=\frac{1}{2}\left ( \frac{M}{M_{Sun}} \right )^{2}}$.

Ignoring the factor $\boldsymbol{ (1/2)^{1/5}=0.87\approx 1}$, we obtain the radius of the sphere of influence:

$\LARGE\boldsymbol{ r_{soi}=R\left ( \frac{M}{M_{Sun}} \right )^{\frac{2}{5}}}$.

(3)

10.  a.  Philae was released from Rosetta when comet 67P was about 3 AU from the Sun, from an initial distance 22.5 km from the comet’s CM.  Show that the spacecraft was well within the comet’s sphere of influence at that time.  (Ans:  $\boldsymbol{ r_{soi}=54\textup{ km}}$)

b.  At perihelion (August 13, 2015), 67P was 1.24 AU from the Sun.  To escape possible damage from the comet’s efflux, Rosetta was removed to a distance of 330 km.  Was it in an unpowered orbit about 67P at that time?  (Ans: No, $\boldsymbol{ r_{soi}=22\textup{ km}}$.  See Ref. 5, fast-forward to August 2015.)

## C. References

1a.  “Modeling Rosetta’s final descent,” Philip Blanco, The Physics Teacher 54, 516 (2016)

2a.  If students are familiar with polar coordinates, substitute $\boldsymbol{ \hat{r}}$ and $\boldsymbol{ \hat{\theta}}$ for $\boldsymbol{ \hat{i}}$ and $\boldsymbol{ \hat{j}}$.