# 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}$.