<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Throwing Axes
Modeling and Simulation in Python
Copyright 2021 Allen Downey
License
Step1:  My favorite event at Lumberjack Competitions is axe throwing.  The axes used for this event typically weigh 1.5 to 2 kg, with handles roughly 0.7 m long.  They are thrown overhead at a target typically 6 m away and 1.5 m off the ground.  Normally, the axe makes one full rotation in the air to hit the target blade first, with the handle close to vertical.
Here's a version of make_system that sets the initial conditions.
The state variables are x, y, theta, vx, vy, omega, where theta is the orientation (angle) of the axe in radians and omega is the angular velocity in radians per second.
I chose initial conditions based on videos of axe throwing.
Step2:  As a simple starting place, I ignore drag, so vx and omega are constant, and ay is just -g.
Step3:  As always, let's test the slope function with the initial conditions.
Step4:  And then run the simulation.
Step5:  Visualizing the results
The simplest way to visualize the results is to plot the state variables as a function of time.
Step6:  We can plot the velocities the same way.
Step7:  Another way to visualize the results is to plot y versus x.  The result is the trajectory through the plane of motion.
Step8:  Animation
Animating this system is a little more complicated, if we want to show the shape and orientation of the axe.
It is useful to construct a frame with $\hat{r}$ along the handle of the axe and $\hat{\theta}$ perpendicular.
Step9:  Now we're ready to animate the results.  The following figure shows the frame and the labeled points A, B, C, and D.
Step10:  During the animation, the parts of the axe seem to slide around relative to each other.  I think that's because the lines and circles get rounded off to the nearest pixel.
Here's the final state of the axe at the point of impact (assuming the target is 8 m away).
<ASSISTANT_TASK:>
Python Code:
# install Pint if necessary
try:
    import pint
except ImportError:
    !pip install pint
# download modsim.py if necessary
from os.path import exists
filename = 'modsim.py'
if not exists(filename):
    from urllib.request import urlretrieve
    url = 'https://raw.githubusercontent.com/AllenDowney/ModSim/main/'
    local, _ = urlretrieve(url+filename, filename)
    print('Downloaded ' + local)
# import functions from modsim
from modsim import *
Explanation: Throwing Axes
Modeling and Simulation in Python
Copyright 2021 Allen Downey
License: Creative Commons Attribution-NonCommercial-ShareAlike 4.0 International
End of explanation
g = 9.8       # m/s**2,
mass = 1.5    # kg,
length = 0.7  # m
x = 0      # m
y = 2      # m
vx = 8     # m / s
vy = 4     # m / s
theta = 2            # rad
omega = -7           # rad / s
init = State(x=x, y=y, 
             vx=vx, vy=vy, 
             theta=theta, omega=omega)
system = System(init=init, t_end=1)
Explanation: My favorite event at Lumberjack Competitions is axe throwing.  The axes used for this event typically weigh 1.5 to 2 kg, with handles roughly 0.7 m long.  They are thrown overhead at a target typically 6 m away and 1.5 m off the ground.  Normally, the axe makes one full rotation in the air to hit the target blade first, with the handle close to vertical.
Here's a version of make_system that sets the initial conditions.
The state variables are x, y, theta, vx, vy, omega, where theta is the orientation (angle) of the axe in radians and omega is the angular velocity in radians per second.
I chose initial conditions based on videos of axe throwing.
End of explanation
def slope_func(t, state, system):
    x, y, vx, vy, theta, omega = state
    ax = 0
    ay = -g
    alpha = 0
    return vx, vy, ax, ay, omega, alpha
Explanation: As a simple starting place, I ignore drag, so vx and omega are constant, and ay is just -g.
End of explanation
slope_func(0, system.init, system)
Explanation: As always, let's test the slope function with the initial conditions.
End of explanation
results, details = run_solve_ivp(system, slope_func)
details
results.tail()
Explanation: And then run the simulation.
End of explanation
def plot_position(results):
    results.x.plot(label='x')
    results.y.plot(label='y')
    decorate(xlabel='Time (s)',
             ylabel='Position (m)')
    
plot_position(results)
Explanation: Visualizing the results
The simplest way to visualize the results is to plot the state variables as a function of time.
End of explanation
def plot_velocity(results):
    results.vx.plot(label='vx')
    results.vy.plot(label='vy')
    decorate(xlabel='Time (s)',
             ylabel='Velocity (m/s)')
    
plot_velocity(results)
results.theta.plot(label='theta', color='C2')
decorate(xlabel='Time (s)',
         ylabel='Angle (radian)')
results.omega.plot(label='omega', color='C2')
decorate(xlabel='Time (s)',
         ylabel='Angular velocity (rad/s)')
Explanation: We can plot the velocities the same way.
End of explanation
def plot_trajectory(results):
    x = results.x
    y = results.y
    make_series(x, y).plot(label='trajectory')
    
    decorate(xlabel='x position (m)',
             ylabel='y position (m)')
    
plot_trajectory(results)
Explanation: Another way to visualize the results is to plot y versus x.  The result is the trajectory through the plane of motion.
End of explanation
def make_frame(theta):
    x, y = pol2cart(theta, 1)
    rhat = Vector(x, y)
    that = vector_perp(rhat)
    return rhat, that
theta = 1
rhat, that = make_frame(theta)
rhat
that
np.dot(rhat, that)
O = Vector(0, 0)
plot_segment(O, rhat)
plot_segment(O, that)
plt.axis('equal');
Explanation: Animation
Animating this system is a little more complicated, if we want to show the shape and orientation of the axe.
It is useful to construct a frame with $\hat{r}$ along the handle of the axe and $\hat{\theta}$ perpendicular.
End of explanation
l1 = 0.6
l2 = 0.1
def draw_func(t, state):   
    x, y, vx, vy, theta, omega = state
    P = Vector(x, y)
    rhat, that = make_frame(theta)
    
    # plot the handle
    A = P - l1 * rhat
    B = P + l2 * rhat
    plot_segment(A, B, color='red')
    # plot the axe head
    C = B + l2 * that
    D = B - l2 * that
    plot_segment(C, D, color='black', linewidth=10)
    # plot the COG
    plt.plot(x, y, 'bo')
    decorate(xlabel='x position (m)',
             ylabel='y position (m)',
             xlim=[-0.3, 8.2],
             ylim=[0, 5]
            )
Explanation: Now we're ready to animate the results.  The following figure shows the frame and the labeled points A, B, C, and D.
End of explanation
draw_func(0, system.init)
animate(results, draw_func)
Explanation: During the animation, the parts of the axe seem to slide around relative to each other.  I think that's because the lines and circles get rounded off to the nearest pixel.
Here's the final state of the axe at the point of impact (assuming the target is 8 m away).
End of explanation
<END_TASK>