Modelling the trajectory of a badminton shuttlecock

%matplotlib widget
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import IPython
from scipy import integrate

v0 = 493.0 / 3.6  # Initial velocity [m/s]
A = 4.0e-3  # Shuttlecock cross area [m**2]
cx = 0.62  # Drag coefficient []
m = 4.0e-2  # Shuttlecock mass [kg]
rho = 1.225  # Air density [kg/m**3]
g = 9.81  # Gravity [m/s**2]

def derivative(X, t):
    """
    Target ODE: Newton's second law
    """
    x, y, vx, vy = X
    v = (vx**2 + vy**2) ** 0.5
    Tx, Ty = vx / v, vy / v
    ax = -0.5 * rho * v**2 * A * cx * Tx / m
    ay = -0.5 * rho * v**2 * A * cx * Ty / m - g
    return np.array([vx, vy, ax, ay])


x0, y0 = 0.0, 0.0
theta0 = 45.0
X0 = [x0, y0, v0 * np.cos(np.radians(theta0)), v0 * np.sin(np.radians(theta0))]
t = np.linspace(0.0, 10.0, 200)
sol = integrate.odeint(derivative, X0, t)
out = pd.DataFrame(sol, columns=["x", "y", "vx", "vy"])
out.head()

plt.figure()
plt.plot(out.x, out.y)
plt.grid()
plt.ylim(0.0, 50.0)
plt.xlabel("Position, $x$")
plt.ylabel("Position, $y$")
plt.show()

thetas = [0.0, 10.0, 15.0, 20.0, 30.0, 45.0, 60.0, 80.0, 85.0]
plt.figure()

for theta0 in thetas:
    x0, y0 = 0.0, 3.0
    X0 = [x0, y0, v0 * np.cos(np.radians(theta0)), v0 * np.sin(np.radians(theta0))]
    t = np.linspace(0.0, 10.0, 1000)
    sol = integrate.odeint(derivative, X0, t)
    out = pd.DataFrame(sol, columns=["x", "y", "vx", "vy"])
    out["t"] = t
    plt.plot(out.x, out.y, label=r"$\theta_0 = $" + "{0}".format(theta0))
plt.legend()
plt.grid()
plt.ylim(0.0, 50.0)
plt.xlabel("Position, $x$")
plt.ylabel("Position, $y$")
plt.show()

%%time
thetas = np.linspace(-180.0, 180.0, 300)
xmax = np.zeros_like(thetas)


for i in range(len(thetas)):
    theta0 = thetas[i]
    x0, y0 = 0.0, 3.0
    X0 = [x0, y0, v0 * np.sin(np.radians(theta0)), -v0 * np.cos(np.radians(theta0))]
    t = np.linspace(0.0, 10.0, 10000)
    sol = integrate.odeint(derivative, X0, t)
    out = pd.DataFrame(sol, columns=["x", "y", "vx", "vy"])
    xmax[i] = out[out.y < 0.0].iloc[0].x

plt.figure()
plt.plot(thetas, xmax)
plt.grid()
plt.xlabel(r"Start angle $\theta_0$")
plt.ylabel(r"Range $x_m$")
plt.show()