A spherical projectile of mass $m$ launched with some initial velocity moves under the influence of two forces: gravity, $\boldsymbol{F}_g = -mg\boldsymbol{\hat{z}}$, and air resistance (drag), $\boldsymbol{F}_D = -\frac{1}{2}c\rho A v^2 \boldsymbol{v}/|\boldsymbol{v}| = -\frac{1}{2}c\rho A v\boldsymbol{v}$, acting in the opposite direction to the projectile's velocity and proportional to the square of that velocity (under most realistic conditions). Here, $c$ is the drag coefficient, $\rho$ the air density, and $A$ the projectile's cross-sectional area.
The relevant equations of motion are therefore:
$$
\begin{align*}
m\ddot{x} &= -k\sqrt{\dot{x}^2 + \dot{z}^2}\dot{x},\\
m\ddot{z} &= -k\sqrt{\dot{x}^2 + \dot{z}^2}\dot{z} - mg,
\end{align*}
$$
where $v = |\boldsymbol{v}| = \sqrt{\dot{x}^2 + \dot{z}^2}$ and $k=\frac{1}{2}c\rho A$. These can be decomposed into the following four first-order ODEs with $u_1 \equiv x, u_2 \equiv \dot{x}, u_3 \equiv z, u_4 \equiv \dot{z}$:
$$
\begin{align*}
\dot{u}_1 &= u_2,\\
\dot{u}_2 &= -\frac{k}{m}\sqrt{u_2^2 + u_4^2}u_2,\\
\dot{u}_3 &= u_4, \\
\dot{u}_4 &= -\frac{k}{m}\sqrt{u_2^2 + u_4^2}u_4 - g.
\end{align*}
$$
The following code integrates this system and identifies two events: hitting the target (the projectile returning to the ground at $z=0$) and reaching its maximum height (at which the z-component of its velocity is zero). We set the additional attribute hit_target.direction = -1 to ensure that hit_target only triggers the event when its return value (the projectile's elevation) goes from positive to negative; otherwise the event would be triggered at launch since $z_0 = 0$. Other possibilities are direction = 1: trigger the event when the return value changes from negative to positive or direction = 0 (the default): the event is triggered when the return value is zero from either direction.
import numpy as np
from scipy.integrate import solve_ivp
import matplotlib.pyplot as plt
# Drag coefficient, projectile radius (m), area (m2) and mass (kg).
c = 0.47
r = 0.05
A = np.pi * r**2
m = 0.2
# Air density (kg.m-3), acceleration due to gravity (m.s-2).
rho_air = 1.28
g = 9.81
# For convenience, define this constant.
k = 0.5 * c * rho_air * A
# Initial speed and launch angle (from the horizontal).
v0 = 50
phi0 = np.radians(65)
def deriv(t, u):
x, xdot, z, zdot = u
speed = np.hypot(xdot, zdot)
xdotdot = -k/m * speed * xdot
zdotdot = -k/m * speed * zdot - g
return xdot, xdotdot, zdot, zdotdot
# Initial conditions: x0, v0_x, z0, v0_z.
u0 = 0, v0 * np.cos(phi0), 0., v0 * np.sin(phi0)
# Integrate up to tf unless we hit the target sooner.
t0, tf = 0, 50
def hit_target(t, u):
# We've hit the target if the z-coordinate is 0.
return u[2]
# Stop the integration when we hit the target.
hit_target.terminal = True
# We must be moving downwards (don't stop before we begin moving upwards!)
hit_target.direction = -1
def max_height(t, u):
# The maximum height is obtained when the z-velocity is zero.
return u[3]
soln = solve_ivp(deriv, (t0, tf), u0, dense_output=True,
events=(hit_target, max_height))
print(soln)
print('Time to target = {:.2f} s'.format(soln.t_events[0][0]))
print('Time to highest point = {:.2f} s'.format(soln.t_events[1][0]))
# A fine grid of time points from 0 until impact time.
t = np.linspace(0, soln.t_events[0][0], 100)
# Retrieve the solution for the time grid and plot the trajectory.
sol = soln.sol(t)
x, z = sol[0], sol[2]
print('Range to target, xmax = {:.2f} m'.format(x[-1]))
print('Maximum height, zmax = {:.2f} m'.format(max(z)))
plt.plot(x, z)
plt.xlabel('x /m')
plt.ylabel('z /m')
plt.show()
The output is:
Time to target = 6.34 s
Time to highest point = 2.79 s
Range to target, xmax = 64.12 m
Maximum height, zmax = 49.42 m
and the plot created is shown in below.
