# Projectile motion
Let's apply our ODE solvers to some problems involving balls and projectiles.

The `integrators.py` file from section 11 is used here (and named `ode.py`).


```python
import numpy as np
import ode
```


```python
%matplotlib inline
import matplotlib.pyplot as plt
```

## Projectile with linear air-resistance
Linear drag force

$$
\mathbf{F}_1 = -b_1 \mathbf{v}
$$

Equations of motion with force due to gravity $\mathbf{g} = -g \hat{\mathbf{e}}_y$

\begin{align}
\frac{d\mathbf{r}}{dt} &= \mathbf{v}\\
\frac{d\mathbf{v}}{dt} &= - g \hat{\mathbf{e}}_y -\frac{b_1}{m} \mathbf{v} 
\end{align}

Bring into standard ODE form for

$$
\frac{d\mathbf{y}}{dt} = \mathbf{f}(t, \mathbf{y})
$$

as

$$
\mathbf{y} = \begin{pmatrix}
x\\
y\\
v_x\\
v_y
\end{pmatrix}, \quad
\mathbf{f} = \begin{pmatrix}
v_x\\
v_y\\
-\frac{b_1}{m} v_x\\
-g -\frac{b_1}{m} v_y
\end{pmatrix}
$$

(Based on Wang 2016, Ch 3.3.1)

**Complete the code**
- Formulate the function `f()` for the standard ODE form
- Set up the integration loop: 
  - only integrate until the particle hits ground, i.e. while $y ≥ 0$.
  - choose an appropriate ODE solver from `ode.py`


```python
def simulate(v0, h=0.01, b1=0.2, g=9.81, m=0.5):
    
    def f(t, y):
        # y = [x, y, vx, vy]
        raise NotImplementedError

    vx, vy = v0
    t = 0
    positions = []
    y = np.array([0, 0, vx, vy], dtype=np.float64)
    
    # add integration loop here
        
    return np.array(positions)

def initial_v(v, theta):
    x = np.deg2rad(theta)
    return v * np.array([np.cos(x), np.sin(x)])
```

Run for an initial velocity of 200 m/s and an angle of $\theta = 30^\circ$:


```python
r = simulate(initial_v(200, 30), h=0.01, b1=1)
```

Plot the motion in the $x$-$y$ plane.


```python
# plot
```

Plot for different starting angles:


```python
for angle in (5, 7.5, 10, 20, 30, 45):
    # ... complete
```

## Simple Baseball physics

- quadratic air resistance (with velocity-dependent drag coefficient)
- Magnus force due to spin


### Quadratic air resistance
Occurs at high Reynolds numbers, i.e., turbulent flow. Only approximate:

$$
\mathbf{F}_2 = -b_2 v \mathbf{v}
$$

### Magnus effect 

**Magnus effect**: airflow is changed around a spinning object. The Magnus force is

$$
\mathbf{F}_M = \alpha \boldsymbol{\omega} \times \mathbf{v}
$$

where $\boldsymbol{\omega}$ is the ball's angular velocity in rad/s (e.g., 200/s for a baseball).

For a sphere the proportionality constant $\alpha$ can be written

$$
\mathbf{F}_M = \frac{1}{2} C_L \rho A \frac{v}{\omega} \boldsymbol{\omega} \times \mathbf{v}
$$

where $C_L$ is the lift coefficient, $\rho$ the air density, $A$ the ball's cross section. (Advantage of defining $C_L$ this way: when spin and velocity are perpendicular, the Magnus force is simply $F_M = \frac{1}{2} C_L \rho A v^2$.)

$C_L$ is mainly a function of the *spin parameter*

$$
S = \frac{r\omega}{v}
$$

with the radius $r$ of the ball. In general we write

$$
\mathbf{F}_M = \frac{1}{2} C_L  \frac{\rho A r}{S} \boldsymbol{\omega} \times \mathbf{v}
$$

For a baseball, experimental data show approximately a power law dependence of $C_L$ on $S$

$$
C_L = 0.62 \times S^{0.7}
$$

All together:

\begin{align}
\mathbf{F}_M &= \alpha\ \boldsymbol{\omega} \times \mathbf{v}\\
v &= \sqrt{\mathbf{v}}\\
S &= \frac{r\omega}{v}\\
C_L &= 0.62 \times S^{0.7}\\
\alpha &= \frac{1}{2} C_L  \frac{\rho A r}{S}
\end{align}


### Equations of motion

\begin{align}
\frac{d\mathbf{r}}{dt} &= \mathbf{v}\\
\frac{d\mathbf{v}}{dt} &= -g \hat{\mathbf{e}}_y \mathbf{v} -\frac{b_2}{m} v \mathbf{v} + \alpha\ \boldsymbol{\omega} \times \mathbf{v}\\
\end{align}

(quadratic drag $-\frac{b_2}{m} v \mathbf{v}$ included.)


### Baseball simulation 

Implement the full baseball equations of motions:
- gravity $a_\text{gravity}$
- quadratic drag $a_\text{drag}$
- Magnus effect $a_\text{Magnus}$

For the cross product you can look at [numpy.cross()](https://docs.scipy.org/doc/numpy/reference/generated/numpy.cross.html).


```python
def C_L(S):
    return 0.62 * S**0.7

def simulate_baseball(v0, omega=200.*np.array([0,1,1]), r0=np.array([0, 2.]),
                      h=0.01, b2=0.0013, g=9.81, rho = 1.225,
                      r=0.07468/2, m=0.14883, R_homeplate=18.4):

    # make sure that omega is a numpy array
    omega = np.asarray(omega)
    
    # all SI units (kg, m)
    # air density rho in kg/m^3
    
    domega = np.linalg.norm(omega)
    A = np.pi*r**2
    rhoArm = rho * A * r / m
    
    # internally, use 3d coordinates [x,y,z];
    # y = [x, y, z, vx, vy, vz]
    
    a_gravity = np.array([0, -g, 0])

    def f(t, y):
        # y = [x, y, z, vx, vy, vz]
        v = y[3:]
        dv = np.linalg.norm(v)
        # COMPLETE
        # 1. acceleration due to drag
        # 2. acceleration due to Magnus effect
        # 3. acceleration due to gravity (a_gravity)
        
        # need to return array f of length 6!
        raise NotImplementedError

    x0, y0 = r0
    vx, vy = v0
    t = 0
    positions = []
    # initialize 3D!
    y = np.array([x0, y0, 0, vx, vy, 0], dtype=np.float64)
    
 
    # IMPLEMENT integration loop
    # - use ode.rk4()
    # - stop when x >= R_homeplate or y < 0.2 (i.e. cannot be caught)
    
    return np.array(positions)

```

#### Simulate throws 

Simulate baseball throw for initial velocity $\mathbf{v} = (30\,\text{m/s}, 0)$.

Plot x vs y and x vs z (to see curving).

Try out different spins; a good value is $\boldsymbol{\omega} = 200\,\text{rad/s} \times (0, 1, 1)$.

Simulate the baseball throw with
- almost no spin: $\omega = 0.01 \times (0, 0, 1)$
- full forward spin: $\omega = 200 \times (0, 0, 1)$
- sideways spin: $\omega = 200 \times (0, 1, 1)$

#### Plotting

Plot the three scenarios in 2D planes: x-y (side view) and x-z (top view).


```python

```

b987322cbe89303d27532c167f8cf1f353ab8193
ipynb
Jupyter Notebook
12_ODE_applications/12-ODE-applications-students.ipynb
nachrisman/PHY494
bac0dd5a7fe6f59f9e2ccaee56ebafcb7d97e2e7
CC-BY-4.0
12_ODE_applications/12-ODE-applications-students.ipynb
nachrisman/PHY494
bac0dd5a7fe6f59f9e2ccaee56ebafcb7d97e2e7
CC-BY-4.0
12_ODE_applications/12-ODE-applications-students.ipynb
nachrisman/PHY494
bac0dd5a7fe6f59f9e2ccaee56ebafcb7d97e2e7
CC-BY-4.0
Qwen/Qwen-72B
1. YES
2. YES