```python
import matplotlib.pyplot as plt
import numpy as np
from scipy.integrate import solve_ivp

from IPython.display import display
import numpy as np
import matplotlib.pyplot as plt
import ipywidgets as widgets
```

# Motion of a spherical projectile with spin and drag

Author: Dr. James Munroe

Date: March 1, 2020

## Projectile motion

According to Newton's Law of Motion an object in motion such a ball moving through the air as projectile only changes its direction and speed (that is, *acceleration*, $a$) in response to the sum of all the forces $\sum F$ acting on the object.  Mathematically we write

$$ a = \frac{\sum F}{m} $$

where $m$ is the mass of the object.

Consider an $x-y$ coordinate system where the $y$-axis the vertical direction. Ground level is given by $y=0$ and $y$ increases going upward.

A projectile at any point in time is described by its position $(x, y)$ and its velocity $(v_x, v_y)$

To predict the motion of an spinning ball with air resistance we have to consider all of the forces involved.

### Gravity

In simplest case, the only force that is important is the force of gravity which is given by 

$$ F_G = m g $$

and points downward towards and the so-called acceleration due to gravity $g$. Close to the surface of the Earth, $g$ is effectively a constant. In this report, we will use a nominal value of $g = 9.81\;\mbox{m/s}^2$.

## Drag

For an object moving through a air, there is air resistance that lead to a drag force

$$ F_D = \frac{1}{2} C_d  \rho  v^2  A $$

where $\rho$ is the density of air, $v$ is the speed of the object relative to the air, and $A$ is the cross-sectional area of the object.  There is also a drag coefficient, $C_d$, which must be experimentally determined.  For a sphere with a rough surface, this value is is approximately $C_d = 0.47$.

The drag force always act in the direction opposite to velocity of the object.

https://www.grc.nasa.gov/www/k-12/airplane/dragsphere.html



## Lift

A rotating object such as a spinning ball will experience a lift force due to the Magnus effect. This lift force is given by

$$ F_L = C_l  \frac{4}{3}\left(4 \pi^2 r^3 s \rho v\right)$$ 

where $r$ is the radius of the oball, $\rho$ is the density of air, $v$ is the speed of the ball, and $s$ is rate of spinning. The lift coefficient $C_l$ is experimental determined value that depends on viscosity and the boundary layer on the ball. As a reference, a baseball has a value of $C_l = 0.15$  

The lift force act perpendicular to direction of motion for the object and perpendicular to the axis of rotation.  For an object with topspin or backspin this will result in a downward or upward force.  If the ball was spinning around vertical axis, such as in a curveball thrown in baseball, this force can act to the left or right. 

https://www.grc.nasa.gov/www/k-12/airplane/balllift.html

## Equations of Motion

For a simple projectile only under the influence of gravity, its motion is given by the equations

$$\begin{align}
\frac{dv_x}{dt} &= 0 \\
\frac{dv_y}{dt} &= \frac{1}{m}\left(-mg\right)
\end{align}
$$ 
These say that the horizontal component of the velocity never changes (since there are no forces in the horizontal) while the vertical component of the velocity has an acceleration given by $-g$.

The kinematic equations learned in high school solve for the position of object in projectile motion but we are going to solve them dynamically so that we can include the effect of additional forces later on.

When we add in drag and lift forces, the equations become

$$\begin{align}
\frac{dv_x}{dt} &= \frac{1}{m}\left( - F_D \frac{v_x}{v} + F_L\frac{v_y}{v} \right) \\
\frac{dv_y}{dt} &= \frac{1}{m}\left( - F_D \frac{v_y}{v} - F_L\frac{v_x}{v} -mg\right)
\end{align}
$$ 

To solve these equations we form a *state vector* contain in both the position and velocity components.

$$ {\bf s} = (x, y, v_x, v_y)$$


```python
def calc_forces(t, s, params):
    
    x, y, vx, vy = s
    
    m = params['m']
    g = params['g']
    Cd = params['Cd']
    Cl = params['Cl']
    ρ = params['ρ']
    r = params['r']
    ω = params['ω']
    
    A = np.pi * r**2
    v = np.sqrt(vx**2 + vy**2)
    
    FG = m*g
    FG_x = 0
    FG_y = - FG
    
    FD = 0.5 * Cd * ρ * A * v**2
    FD_x = -FD*vx/v
    FD_y = -FD*vy/v
    
    FL = 4/3*Cl*(4*np.pi*r**3 * ω * ρ * v)
    FL_x = FL/v*vy
    FL_y = -FL/v*vx

    
    return { 'FG_x': FG_x, 'FG_y': FG_y,
             'FD_x': FD_x, 'FD_y': FD_y,
             'FL_x': FL_x, 'FL_y': FL_y,
           }

    dxdt = vx
    dydt = vy
    dvxdt = Fnet_x / m
    dvydt = Fnet_y / m
    
    dsdt = [dxdt, dydt, dvxdt, dvydt]
    return dsdt

def projectile_motion(t, s, params):
    
    x, y, vx, vy = s
    
    m = params['m']
    
    forces = calc_forces(t, s, params)
    FG_x = forces['FG_x']
    FG_y = forces['FG_y']
    FD_x = forces['FD_x']
    FD_y = forces['FD_y']
    FL_x = forces['FL_x']
    FL_y = forces['FL_y']
    
    Fnet_x = FG_x + FD_x + FL_x
    Fnet_y = FG_y + FD_y + FL_y
    
    dxdt = vx
    dydt = vy
    dvxdt = Fnet_x / m
    dvydt = Fnet_y / m
    
    dsdt = [dxdt, dydt, dvxdt, dvydt]
    return dsdt


# events to track
def apex(t, s, params): 
    return s[3]
def hit_ground(t, s, params): 
    return s[1]
hit_ground.terminal = True
hit_ground.direction = -1

def plot_solution(t_=0, 
                  x0=0, y0=0, 
                  v0 = 0, # initial speed, m/s
                  θ0 = 0, # initial angle, degrees
                  S = 0, # spin, RPM
                  Cd = 0.5, # drag coefficient
                  Cl = 0.15, # lift coefficient
                 # r = 0.05, # radius, m
                  m = 0.142, # mass, kg
                  show_forces=False,
                  zoom_out=False,
                  
                  **kwargs ):

    r = 7 * 2.54 / 2 / 100 # m
    g = 9.81 # m/s^2
    ρ = 1.2 # density of air
    ω = -S/60*2*np.pi
    tmax = 10
    vx0 = v0 * np.cos(np.deg2rad(θ0))
    vy0 = v0 * np.sin(np.deg2rad(θ0))
    
    params = {'g': g, 'm': m,
              'r' : r,'ρ' : ρ, 'ω': ω,
             'Cd' : Cd, 'Cl': Cl}

    # initial state vector
    s0 = (x0, y0, vx0, vy0)

    # solve equations
    sol = solve_ivp(projectile_motion, 
                [0, tmax], # tspan 
                s0, 
                args=(params,),
                events=(apex, hit_ground), 
                dense_output=True)
  
    try:
        t_apex = sol.t_events[0][0]
        t_ground = sol.t_events[1][0]
    except:
        t_apex = None
        t_ground = tmax

    t = np.linspace(0, t_ground, 100)
    x, y, vx, vy = sol.sol(t)
    
    fig, ax = plt.subplots(figsize=(10, 8))
    plt.plot(x, y, label='Predicted Path')
    
    if t_ > t_ground:
        t_ = t_ground
        
    s_ = sol.sol(t_)
    x_, y_, vx_, vy_ = s_
    plt.plot(x_, y_, 'o')
    
    if show_forces:
        forces = calc_forces(t_, s_, params)
        FG_x = forces['FG_x']
        FG_y = forces['FG_y']
        FD_x = forces['FD_x']
        FD_y = forces['FD_y']
        FL_x = forces['FL_x']
        FL_y = forces['FL_y']
        plt.arrow(x_, y_, FG_x, FG_y, width=0.05, color='red', alpha=0.3)
        plt.annotate('$F_G$', (x_+FG_x, y_+FG_y))
        plt.arrow(x_, y_, FD_x, FD_y, width=0.05, color='blue', alpha=0.3)
        plt.annotate('$F_D$', (x_+FD_x, y_+FD_y))
        if abs(S) > 0:
            plt.arrow(x_, y_, FL_x, FL_y, width=0.05, color='green', alpha=0.3)
            plt.annotate('$F_L$', (x_+FL_x, y_+FL_y))
    
    # solve equations without lift or drag
    params['Cd'] = 0
    params['Cl'] = 0
    sol = solve_ivp(projectile_motion, 
                [0, tmax], # tspan 
                s0, 
                args=(params,),
                events=(apex, hit_ground), 
                dense_output=True)
    try:
        t_apex = sol.t_events[0][0]
        t_ground = sol.t_events[1][0]
    except:
        t_apex = None
        t_ground = tmax
    
    t = np.linspace(0, t_ground, 100)
    x, y, vx, vy = sol.sol(t)
    plt.plot(x, y, '--', alpha=0.6, label='Idealized Parabola')
    
    plt.xlabel('x (m)')
    plt.ylabel('y (m)')
    plt.axis('equal')
    
    plt.axhline(color='k')
    rect=plt.Rectangle((x0-0.75,0), 0.75, 0.8)
    ax.add_patch(rect)
    
    W = 29.25 * 2.54 / 100
    rect=plt.Rectangle((10,0), W, 4, color='red', alpha=0.6)
    ax.add_patch(rect)
    VF_1 = 83.25 * 2.54 / 100
    VF_2 = 113.25 * 2.54 / 100
    VB = 98.25 * 2.54 / 100
    VB_d = 13 / 2 * 2.54 / 100
    
    rect=plt.Rectangle((10, VF_1), W, (VF_2 - VF_1), color='purple',  alpha=0.6)
    ax.add_patch(rect)
    rect=plt.Rectangle((10+0.75*W, VB-VB_d/2), 0.25*W, VB_d, color='yellow',  alpha=0.6)
    ax.add_patch(rect)
    
    if not zoom_out:
        plt.xlim(-2, 12)
        plt.ylim(0, 4)
    
    plt.legend()
    plt.show()
    
w = dict(t_ = widgets.FloatSlider(0, min=0.0, max=2, step=0.01, description='t'),
         x0 = widgets.FloatText(value=0, step=0.1, description='$x_0\;(\mbox{m})$:'),
         y0 = widgets.BoundedFloatText(value=0.99, min=0.0, step=0.01, description='$y_0\;(\mbox{m})$:'),
         θ0 = widgets.BoundedFloatText(value=15, min=0, max=90, step=1, description='$θ_{0}\;(\mbox{degrees})$:'),
         v0 = widgets.FloatText(value=12, step=0.1, description='$v\;(\mbox{m/s})$:'),
         m = widgets.FloatText(value=0.142 , step=0.001, description='$m\;(\mbox{kg})$:'),
         Cd = widgets.BoundedFloatText(value=0.47 , min = 0.0, max = 2.0, step=0.01, description='$C_d$:'),
         S = widgets.FloatText(value=600 ,  step=10, description='ball RPM:'),
         Cl = widgets.BoundedFloatText(value=0.15 , min = 0.0, max = 1.0, step=0.01, description='$C_l$:'),
         show_forces = widgets.Checkbox(value=False, description='Show force vectors'),
         zoom_out = widgets.Checkbox(value=False, description='Zoom out')
                              )

instructions = widgets.Label('Notes: ball spins counter-clockwise')
output = widgets.interactive_output(plot_solution, w)
box = widgets.HBox([output, widgets.VBox([*w.values()])])
app = widgets.VBox([widgets.Label('Projectile Motion of a Spinning Ball with Air Drag'), box, instructions])
display(app)
```


    VBox(children=(Label(value='Projectile Motion of a Spinning Ball with Air Drag'), HBox(children=(Output(), VBo…



```python

```

b8de1ca36572acdb116a1217eb44546310dc11ed
ipynb
Jupyter Notebook
index.ipynb
jmunroe/projectile-motion
c53b7829bac0921e4f8895de1db7370613aeaf15
MIT
index.ipynb
jmunroe/projectile-motion
c53b7829bac0921e4f8895de1db7370613aeaf15
MIT
index.ipynb
jmunroe/projectile-motion
c53b7829bac0921e4f8895de1db7370613aeaf15
MIT
Qwen/Qwen-72B
1. YES
2. YES