Projectile motion

Since there is only acceleration in the vertical direction, the velocity in the horizontal direction is constant, being equal to ${\displaystyle \mathbf {v} _{0}\cos \theta }$. The vertical motion of the projectile is the motion of a particle during its free fall. Here the acceleration is constant, being equal to g.[note 1] The components of the acceleration are:

${\displaystyle a_{x}=0}$,

${\displaystyle a_{y}=-g}$.

Let the projectile be launched with an initial velocity ${\displaystyle \mathbf {v} (0)\equiv \mathbf {v} _{0}}$, which can be expressed as the sum of horizontal and vertical components as follows:

${\displaystyle \mathbf {v} _{0}=v_{0x}\mathbf {\hat {x}} +v_{0y}\mathbf {\hat {y}} }$.

The components ${\displaystyle v_{0x}}$ and ${\displaystyle v_{0y}}$ can be found if the initial launch angle, ${\displaystyle \theta }$, is known:

${\displaystyle v_{0x}=v_{0}\cos(\theta )}$,

${\displaystyle v_{0y}=v_{0}\sin(\theta )}$

The horizontal component of the velocity of the object remains unchanged throughout the motion. The vertical component of the velocity changes linearly,[note 2] because the acceleration due to gravity is constant. The accelerations in the x and y directions can be integrated to solve for the components of velocity at any time t, as follows:

${\displaystyle v_{x}=v_{0}\cos(\theta )}$,

${\displaystyle v_{y}=v_{0}\sin(\theta )-gt}$.

The magnitude of the velocity (under the Pythagorean theorem, also known as the triangle law):

${\displaystyle v={\sqrt {v_{x}^{2}+v_{y}^{2}}}}$.

Displacement and coordinates of parabolic throwing

At any time ${\displaystyle t}$, the projectile's horizontal and vertical displacement are:

${\displaystyle x=v_{0}t\cos(\theta )}$,

${\displaystyle y=v_{0}t\sin(\theta )-{\frac {1}{2}}gt^{2}}$.

The magnitude of the displacement is:

${\displaystyle \Delta r={\sqrt {x^{2}+y^{2}}}}$.

Consider the equations,

${\displaystyle x=v_{0}t\cos(\theta ),y=v_{0}t\sin(\theta )-{\frac {1}{2}}gt^{2}}$.

If t is eliminated between these two equations the following equation is obtained:

${\displaystyle y=\tan(\theta )\cdot x-{\frac {g}{2v_{0}^{2}\cos ^{2}\theta }}\cdot x^{2}}$.

Since g, θ, and v₀ are constants, the above equation is of the form

${\displaystyle y=ax+bx^{2}}$,

in which a and b are constants. This is the equation of a parabola, so the path is parabolic. The axis of the parabola is vertical.

If the projectile's position (x,y) and launch angle (θ or α) are known, the initial velocity can be found solving for v₀ in the aforementioned parabolic equation:

${\displaystyle v_{0}={\sqrt {{x^{2}g} \over {x\sin 2\theta -2y\cos ^{2}\theta }}}}$.

The parabolic trajectory of a projectile can also be expressed in polar coordinates instead of Cartesian coordinates. In this case, the position has the general formula

${\displaystyle r(\phi )={\frac {v_{0}^{2}\cos ^{2}\theta }{g}}\left(-\tan \phi \sec \phi +{\sqrt {(\tan ^{2}\phi +\tan ^{2}\theta )\sec ^{2}\phi }}\right)}$.

In this equation, the origin is the midpoint of the horizontal range of the projectile, and if the ground is flat, the parabolic arc is plotted in the range ${\displaystyle 0\leq \phi \leq \pi }$. This expression can be obtained by transforming the Cartesian equation as stated above by ${\displaystyle y=r\sin \phi }$ and ${\displaystyle x=r\cos \phi }$.

The total time t for which the projectile remains in the air is called the time of flight.

${\displaystyle y=v_{0}t\sin(\theta )-{\frac {1}{2}}gt^{2}}$

After the flight, the projectile returns to the horizontal axis (x-axis), so ${\displaystyle y=0}$.

${\displaystyle 0=v_{0}t\sin(\theta )-{\frac {1}{2}}gt^{2}}$

${\displaystyle v_{0}t\sin(\theta )={\frac {1}{2}}gt^{2}}$

${\displaystyle v_{0}\sin(\theta )={\frac {1}{2}}gt}$

${\displaystyle t={\frac {2v_{0}\sin(\theta )}{g}}}$

Note that we have neglected air resistance on the projectile.

If the starting point is at height y₀ with respect to the point of impact, the time of flight is:

${\displaystyle t={\frac {d}{v\cos \theta }}={\frac {v\sin \theta +{\sqrt {(v\sin \theta )^{2}+2gy_{0}}}}{g}}}$

As above, this expression can be reduced to

${\displaystyle t={\frac {v\sin {\theta }+{\sqrt {(v\sin {\theta })^{2}}}}{g}}={\frac {v\sin {\theta }+v\sin {\theta }}{g}}={\frac {2v\sin {\theta }}{g}}={\frac {2v\sin {(45)}}{g}}={\frac {2v{\frac {\sqrt {2}}{2}}}{g}}={\frac {{\sqrt {2}}v}{g}}}$

if θ is 45° and y₀ is 0.

Maximum height of projectile

The greatest height that the object will reach is known as the peak of the object's motion. The increase in height will last until ${\displaystyle v_{y}=0}$, that is,

${\displaystyle 0=v_{0}\sin(\theta )-gt_{h}}$.

Time to reach the maximum height(h):

${\displaystyle t_{h}={\frac {v_{0}\sin(\theta )}{g}}}$.

For the vertical displacement of the maximum height of projectile:

${\displaystyle h=v_{0}t_{h}\sin(\theta )-{\frac {1}{2}}gt_{h}^{2}}$

${\displaystyle h={\frac {v_{0}^{2}\sin ^{2}(\theta )}{2g}}}$

The maximum reachable height is obtained for θ=90°:

${\displaystyle h_{\mathrm {max} }={\frac {v_{0}^{2}}{2g}}}$

The relation between the range d on the horizontal plane and the maximum height h reached at ${\displaystyle {\frac {t_{d}}{2}}}$ is:

${\displaystyle h={\frac {d\tan \theta }{4}}}$

The maximum distance of projectile

The range and the maximum height of the projectile does not depend upon its mass. Hence range and maximum height are equal for all bodies that are thrown with the same velocity and direction. The horizontal range d of the projectile is the horizontal distance it has traveled when it returns to its initial height (${\displaystyle y=0}$).

${\displaystyle 0=v_{0}t_{d}\sin(\theta )-{\frac {1}{2}}gt_{d}^{2}}$.

Time to reach ground:

${\displaystyle t_{d}={\frac {2v_{0}\sin(\theta )}{g}}}$.

From the horizontal displacement the maximum distance of projectile:

${\displaystyle d=v_{0}t_{d}\cos(\theta )}$,

so[note 3]

${\displaystyle d={\frac {v_{0}^{2}}{g}}\sin(2\theta )}$.

Note that d has its maximum value when

${\displaystyle \sin 2\theta =1}$,

which necessarily corresponds to

${\displaystyle 2\theta =90^{\circ }}$,

or

${\displaystyle \theta =45^{\circ }}$.

Trajectories of projectiles launched at different elevation angles but the same speed of 10 m/s in a vacuum and uniform downward gravity field of 10 m/s². Points are at 0.05 s intervals and length of their tails is linearly proportional to their speed. t = time from launch, T = time of flight, R = range and H = highest point of trajectory (indicated with arrows).

The total horizontal distance (d) traveled.

${\displaystyle d={\frac {v\cos \theta }{g}}\left(v\sin \theta +{\sqrt {(v\sin \theta )^{2}+2gy_{0}}}\right)}$

When the surface is flat (initial height of the object is zero), the distance traveled:[3]

${\displaystyle d={\frac {v^{2}\sin(2\theta )}{g}}}$

Thus the maximum distance is obtained if θ is 45 degrees. This distance is:

${\displaystyle d_{\mathrm {max} }={\frac {v^{2}}{g}}}$

According to the work-energy theorem the vertical component of velocity is:

${\displaystyle v_{y}^{2}=(v_{0}\sin \theta )^{2}-2gy}$.

These formulae ignore aerodynamic drag and also assume that the landing area is at uniform height 0.

The "angle of reach" is the angle (θ) at which a projectile must be launched in order to go a distance d, given the initial velocity v.

${\displaystyle \sin(2\theta )={\frac {gd}{v^{2}}}}$

There are two solutions:

${\displaystyle \theta ={\frac {1}{2}}\arcsin \left({\frac {gd}{v^{2}}}\right)}$ (shallow trajectory)

and

${\displaystyle \theta ={\frac {1}{2}}\arccos \left({\frac {gd}{v^{2}}}\right)}$ (steep trajectory)

Vacuum trajectory of a projectile for different launch angles. Launch speed is the same for all angles, 50 m/s if "g" is 10 m/s².

To hit a target at range x and altitude y when fired from (0,0) and with initial speed v the required angle(s) of launch θ are:

${\displaystyle \tan \theta ={\left({\frac {v^{2}\pm {\sqrt {v^{4}-g(gx^{2}+2yv^{2})}}}{gx}}\right)}}$

The two roots of the equation correspond to the two possible launch angles, so long as they aren't imaginary, in which case the initial speed is not great enough to reach the point (x,y) selected. This formula allows one to find the angle of launch needed without the restriction of ${\displaystyle y=0}$.

One can also ask what launch angle allows the lowest possible launch velocity. This occurs when the two solutions above are equal, implying that the quantity under the square root sign is zero. This requires solving a quadratic equation for ${\displaystyle v^{2}}$, and we find

${\displaystyle v^{2}/g=y+{\sqrt {y^{2}+x^{2}}}.}$

This gives

${\displaystyle \theta =\arctan \left(y/x+{\sqrt {y^{2}/x^{2}+1}}\right).}$

If we denote the angle whose tangent is y/x by α, then

${\displaystyle \tan \theta ={\frac {\sin \alpha +1}{\cos \alpha }}}$

${\displaystyle \tan(\pi /2-\theta )={\frac {\cos \alpha }{\sin \alpha +1}}}$

${\displaystyle \cos ^{2}(\pi /2-\theta )={\frac {1}{2}}(\sin \alpha +1)}$

${\displaystyle 2\cos ^{2}(\pi /2-\theta )-1=\cos(\pi /2-\alpha )}$

This implies

${\displaystyle \theta =\pi /2-{\frac {1}{2}}(\pi /2-\alpha ).}$

In other words, the launch should be at the angle halfway between the target and Zenith (vector opposite to Gravity)

The length of the parabolic arc traced by a projectile L, given that the height of launch and landing is the same and that there is no air resistance, is given by the formula:

${\displaystyle L={\frac {v_{0}^{2}}{2g}}\left(2\sin \theta +\cos ^{2}\theta \cdot \ln {\frac {1+\sin \theta }{1-\sin \theta }}\right)={\frac {v_{0}^{2}}{g}}\left(\sin \theta +\cos ^{2}\theta \cdot \tanh ^{-1}(\sin \theta )\right)}$

where ${\displaystyle v_{0}}$ is the initial velocity, ${\displaystyle \theta }$ is the launch angle and ${\displaystyle g}$ is the acceleration due to gravity as a positive value. The expression can be obtained by evaluating the arc length integral for the height-distance parabola between the bounds initial and final displacements (i.e. between 0 and the horizontal range of the projectile) such that:

${\displaystyle L=\int _{0}^{\mathrm {range} }{\sqrt {1+\left({\frac {\mathrm {d} y}{\mathrm {d} x}}\right)^{2}}}\,\mathrm {d} x=\int _{0}^{v_{0}^{2}\sin(2\theta )/g}{\sqrt {1+\left(-{\frac {g}{v_{0}^{2}\cos ^{2}\theta }}x+\tan \theta \right)^{2}}}\,\mathrm {d} x}$.

Stokes drag, where ${\displaystyle \mathbf {F_{air}} \propto \mathbf {v} }$, only applies at very low speed in air, and is thus not the typical case for projectiles. However, the linear dependence of ${\displaystyle F_{\mathrm {air} }}$ on ${\displaystyle v}$ causes a very simple differential equation of motion

${\displaystyle {\frac {\mathrm {d} }{\mathrm {d} t}}{\begin{pmatrix}v_{x}\\v_{y}\end{pmatrix}}={\begin{pmatrix}-\mu \,v_{x}\\-g-\mu \,v_{y}\end{pmatrix}}}$

in which the two cartesian components become completely independent, and thus easier to solve.[5] Here, ${\displaystyle v_{0}}$,${\displaystyle v_{x}}$ and ${\displaystyle v_{y}}$ will be used to denote the initial velocity, the velocity along the direction of x and the velocity along the direction of y, respectively. The mass of the projectile will be denoted by m, and ${\displaystyle \mu$ :=k/m} . For the derivation only the case where ${\displaystyle 0^{o}\leq \theta \leq 180^{o}}$ is considered. Again, the projectile is fired from the origin (0,0).

${\displaystyle x(t)={\frac {v_{x0}}{\mu }}\left(1-e^{-\mu t}\right)}$ (1b)

${\displaystyle y(t)=-{\frac {g}{\mu }}t+{\frac {1}{\mu }}\left(v_{y0}+{\frac {g}{\mu }}\right)\left(1-e^{-\mu t}\right)}$ (3b)

${\displaystyle t={\frac {1}{\mu }}\left(1+{\frac {\mu }{g}}v_{y0}+W\left(-\left(1+{\frac {\mu }{g}}v_{y0}\right)e^{-\left(1+{\frac {\mu }{g}}v_{y0}\right)}\right)\right)}$.

Trajectories of a skydiver in air with Newton drag

The most typical case of air resistance, for the case of Reynolds numbers above about 1000 is Newton drag with a drag force proportional to the speed squared, ${\displaystyle F_{\mathrm {air} }=-kv^{2}}$. In air, which has a kinematic viscosity around ${\displaystyle 0.15\,\mathrm {cm^{2}/s} }$, this means that the product of speed and diameter must be more than about ${\displaystyle 0.015\,\mathrm {m^{2}/s} }$.

Unfortunately, the equations of motion can not be easily solved analytically for this case. Therefore, a numerical solution will be examined.

The following assumptions are made:

• Constant gravitational acceleration
• Air resistance is given by the following drag formula,

${\displaystyle \mathbf {F_{D}} =-{\tfrac {1}{2}}c\rho A\,v\,\mathbf {v} }$

Where:

• F_D is the drag force
• c is the drag coefficient
• ρ is the air density
• A is the cross sectional area of the projectile
• μ = k/m = cρA/(2m)

Even though the general case of a projectile with Newton drag cannot be solved analytically, some special cases can. Here we denote the terminal velocity in free-fall as ${\displaystyle v_{\infty }={\sqrt {g/\mu }}}$ and the characteristic settling time constant ${\displaystyle t_{f}=1/{\sqrt {g\mu }}}$.

• Near-horizontal motion: In case the motion is almost horizontal, ${\displaystyle |v_{x}|\gg |v_{y}|}$, such as a flying bullet, the vertical velocity component has very little influence on the horizontal motion. In this case:

${\displaystyle {\dot {v}}_{x}(t)=-\mu \,v_{x}^{2}(t)}$

${\displaystyle v_{x}(t)={\frac {1}{1/v_{x,0}+\mu \,t}}}$

${\displaystyle x(t)={\frac {1}{\mu }}\ln(1+\mu \,v_{x,0}\cdot t)}$

The same pattern applies for motion with friction along a line in any direction, when gravity is negligible. It also applies when vertical motion is prevented, such as for a moving car with its engine off.

• Vertical motion upward:

${\displaystyle {\dot {v}}_{y}(t)=-g-\mu \,v_{y}^{2}(t)}$

${\displaystyle v_{y}(t)=v_{\infty }\tan {\frac {t_{\mathrm {peak} }-t}{t_{f}}}}$

${\displaystyle y(t)=y_{\mathrm {peak} }+{\frac {1}{\mu }}\ln \left(\cos {\frac {t_{\mathrm {peak} }-t}{t_{f}}}\right)}$

A projectile can not rise longer than ${\displaystyle t_{\mathrm {rise} }={\frac {\pi }{2}}t_{f}}$ vertically before it reaches the peak.

• Vertical motion downward:[7]

${\displaystyle {\dot {v}}_{y}(t)=-g+\mu \,v_{y}^{2}(t)}$

${\displaystyle v_{y}(t)=-v_{\infty }\tanh {\frac {t-t_{\mathrm {peak} }}{t_{f}}}}$

${\displaystyle y(t)=y_{\mathrm {peak} }-{\frac {1}{\mu }}\ln \left(\cosh {\frac {t-t_{\mathrm {peak} }}{t_{f}}}\right)}$

After a time ${\displaystyle t_{f}}$, the projectile reaches almost terminal velocity ${\displaystyle -v_{\infty }}$.

The approach will be to formulate integral expressions that can be evaluated numerically. All variables will then be expressed in terms of a parameter ${\displaystyle Q}$.

${\displaystyle t(Q)=-{\frac {1}{\sqrt {g\mu }}}\int _{Q_{0}}^{Q}{\frac {\cosh {\tilde {Q}}}{\sqrt {\lambda }}}\,\mathrm {d} {\tilde {Q}}}$

(I)

Also; ${\displaystyle {\frac {\mathrm {d} x}{\mathrm {d} Q}}={\frac {{\Bigl (}{\frac {\mathrm {d} x}{\mathrm {d} t}}{\Bigr )}}{{\Bigl (}{\frac {\mathrm {d} Q}{\mathrm {d} t}}{\Bigr )}}}=-{\frac {1}{\mu }}{\frac {\cosh {Q}}{\lambda }}}$

${\displaystyle x(Q)=x_{0}-{\frac {1}{\mu }}\int _{Q_{0}}^{Q}{\frac {\cosh {\tilde {Q}}}{\lambda }}\,\mathrm {d} {\tilde {Q}}}$

(J)

And; ${\displaystyle {\frac {\mathrm {d} y}{\mathrm {d} Q}}={\frac {{\Bigl (}{\frac {\mathrm {d} y}{\mathrm {d} t}}{\Bigr )}}{{\Bigl (}{\frac {\mathrm {d} Q}{\mathrm {d} t}}{\Bigr )}}}=-{\frac {1}{2\mu }}{\frac {\sinh {2Q}}{\lambda }}}$

${\displaystyle y(Q)=y_{0}-{\frac {1}{2\mu }}\int _{Q_{0}}^{Q}{\frac {\sinh {2{\tilde {Q}}}}{\lambda }}\,\mathrm {d} {\tilde {Q}}}$

(K)

Determine the time of flight ${\displaystyle T}$ by setting ${\displaystyle y}$ to ${\displaystyle 0}$ in equation (K) above. Solve for the value of the variable ${\displaystyle Q}$.

${\displaystyle y_{0}-{\frac {1}{2\mu }}\int _{Q_{0}}^{Q_{T}}{\frac {\sinh {2{\tilde {Q}}}}{\lambda }}\,\mathrm {d} {\tilde {Q}}=0}$

(L)

Equation (I) with ${\displaystyle Q_{T}}$ substituted for ${\displaystyle Q}$ gives;

${\displaystyle T=-{\frac {1}{\sqrt {g\mu }}}\int _{Q_{0}}^{Q_{T}}{\frac {\cosh {\tilde {Q}}}{\sqrt {\lambda }}}\,\mathrm {d} {\tilde {Q}}}$

(M)

Equation (J) gives the horizontal range ${\displaystyle R}$ as;

${\displaystyle R=x_{0}-{\frac {1}{\mu }}\int _{Q_{0}}^{Q_{T}}{\frac {\cosh {\tilde {Q}}}{\lambda }}\,\mathrm {d} {\tilde {Q}}}$

(N)

At the highest point of the projectile path ${\displaystyle \Psi =0}$, and ${\displaystyle Q=0}$, giving the maximum height ${\displaystyle H}$ from equation (K) as;

${\displaystyle H=y_{0}-{\frac {1}{2\mu }}\int _{Q_{0}}^{0}{\frac {\sinh {2{\tilde {Q}}}}{\lambda }}\,\mathrm {d} {\tilde {Q}}}$

(O)

A projectile motion with drag can be computed generically by numerical integration of the ordinary differential equation, for instance by applying a reduction to a first-order system. The equation to be solved is

${\displaystyle {\frac {\mathrm {d} }{\mathrm {d} t}}{\begin{pmatrix}x\\y\\v_{x}\\v_{y}\end{pmatrix}}={\begin{pmatrix}v_{x}\\v_{y}\\-\mu \,v_{x}{\sqrt {v_{x}^{2}+v_{y}^{2}}}\\-g-\mu \,v_{y}{\sqrt {v_{x}^{2}+v_{y}^{2}}}\end{pmatrix}}}$.

The following computer program in the form of a Python script demonstrates such a simulation, where the projectile is modelled as a baseball (parameters from [8]). The script uses the libraries numpy (for arrays), scipy (for numerical integration of the ordinary differential equation, and for root-finding by Newton's method) and matplotlib (for plotting).

#!/usr/bin/env python3
from math import *
import numpy as np
from scipy.integrate import odeint
from scipy.optimize import newton
import matplotlib.pyplot as plt

def projectile_motion(g, mu, xy0, vxy0, tt):
    # use a four-dimensional vector function vec = [x, y, vx, vy]
    def dif(vec, t):
        # time derivative of the whole vector vec
        v = sqrt(vec[2] ** 2 + vec[3] ** 2)
        return [vec[2], vec[3], -mu * v * vec[2], -g - mu * v * vec[3]]

    # solve the differential equation numerically
    vec = odeint(dif, [xy0[0], xy0[1], vxy0[0], vxy0[1]], tt)
    return vec[:, 0], vec[:, 1], vec[:, 2], vec[:, 3]  # return x, y, vx, vy

# Parameters of projectile (modelled after a baseball)
g       = 9.81         # Acceleration due to gravity (m/s^2)
rho_air = 1.29         # Air density (kg/m^3)
v0      = 44.7         # Initial velocity (m/s)
alpha0  = radians(75)  # Launch angle (deg.)
m       = 0.145        # Mass of projectile (kg)
cD      = 0.5          # Drag coefficient (spherical projectile)
r       = 0.0366       # Radius of projectile (m)
mu = 0.5 * cD * (pi * r ** 2) * rho_air / m

# Initial position and launch velocity
x0, y0 = 0.0, 0.0
vx0, vy0 = v0 * cos(alpha0), v0 * sin(alpha0)

T_peak = newton(lambda t: projectile_motion(g, mu, (x0, y0), (vx0, vy0), [0, t])[3][1], 0)
y_peak = projectile_motion(g, mu, (x0, y0), (vx0, vy0), [0, T_peak])[1][1]
T = newton(lambda t: projectile_motion(g, mu, (x0, y0), (vx0, vy0), [0, t])[1][1], 2 * T_peak)
t = np.linspace(0, T, 501)
x, y, vx, vy = projectile_motion(g, mu, (x0, y0), (vx0, vy0), t)

print("Time of flight: {:.1f} s".format(T))        # returns  6.6 s
print("Horizontal range: {:.1f} m".format(x[-1]))  # returns 43.7 m
print("Maximum height: {:.1f} m".format(y_peak))   # returns 53.4 m

# Plot of trajectory
fig, ax = plt.subplots()
ax.plot(x, y, "r-", label="Numerical")
ax.set_title(r"Projectile path")
ax.set_aspect("equal")
ax.grid(b=True)
ax.legend()
ax.set_xlabel("$x$ (m)")
ax.set_ylabel("$y$ (m)")
plt.savefig("01 Path.png")

fig, ax = plt.subplots()
ax.plot(t, vx, "b-", label="$v_x$")
ax.set_title(r"Horizontal velocity component")
ax.grid(b=True)
ax.legend()
ax.set_xlabel("$t$ (s)")
ax.set_ylabel("$v_x$ (m/s)")
plt.savefig("02 Horiz vel.png")

fig, ax = plt.subplots()
ax.plot(t, vy, "b-", label="$v_y$")
ax.set_title(r"Vertical velocity component")
ax.grid(b=True)
ax.legend()
ax.set_xlabel("$t$ (s)")
ax.set_ylabel("$v_y$ (m/s)")
plt.savefig("03 Vert vel.png")