# /usr/bin/env python
import matplotlib.pyplot as plt
import matplotlib.animation as ani
from math import *
[docs]def draw_trajectory(ball):
"""
Draws the trajectory of the ball as a line.
Args:
ball (Ball): the ball whose trajectory will be drawn
"""
# set limits for the image
ax = plt.axes()
ax.set_aspect('equal')
ax.set(xlim=[-2,10], ylim=[-2,10])
# plot horizontal line to denote ground
ax.plot([-2,10],[0,0],'b-')
# plot trajectory
ax.plot( ball.x_trajectory, ball.y_trajectory, 'r-' )
plt.show()
[docs]def draw_motion_blur_image(ball):
"""
Draws the trajectory of the ball as a series or points.
Args:
ball (Ball): the ball whose trajectory will be drawn
"""
# set limits for the image
ax = plt.axes()
ax.set_aspect('equal')
ax.set(xlim=[-2,10], ylim=[-2,10])
# plot horizontal line to denote ground
ax.plot([-2,10],[0,0],'b-')
# plot the ball as a circle
ax.plot( ball.x_trajectory, ball.y_trajectory, 'ro', markersize=6, alpha=0.3 )
plt.show()
[docs]def draw_animation_frame(frame, xs, ys):
"""
Draws the system at a certain time. Used for animation.
Args:
frame (int): index of the time step to draw
xs (list): array of ball x-coordinates
ys (list): array of ball y-coordinates
"""
# clear the image
plt.clf()
# set limits for the image
ax = plt.axes()
ax.set_aspect('equal')
ax.set(xlim=[-2,10], ylim=[-2,10])
# plot horizontal line to denote ground
ax.plot([-2,10],[0,0],'b-')
# plot the ball as a circle
ax.plot( xs[frame], ys[frame], 'ro', markersize=6 )
[docs]def animate(ball):
"""
Animate the motion of the ball.
Args:
ball (Ball): the ball whose trajectory will be drawn
"""
nframes = len(ball.x_trajectory)
print("animating "+str(nframes)+" frames")
fig = plt.figure()
motion = ani.FuncAnimation(fig, draw_animation_frame, nframes, interval=10,
fargs=(ball.x_trajectory, ball.y_trajectory) )
plt.show()
[docs]class Ball:
"""
A center-of-mass model for a ball.
The model does not include rotation of the ball.
Args:
x (float): initial x-coordinate
y (float): initial y-coordinate
vx (float): initial velocity x-component
vy (float): initial velocity y-component
mass (float): ball mass
radius (float): ball radius, needed for bouncing
"""
def __init__(self, x, y, vx, vy, mass = 1.0, radius = 0.1):
# position of the ball
self.x = x
self.y = y
# velocity of the ball
self.vx = vx
self.vy = vy
# force acting on the ball
self.fx = 0.0
self.fy = 0.0
# mass of the ball
self.mass = mass
# radius of the ball
self.radius = radius
# lists for storing the coordinates of the ball as time advances
self.x_trajectory = [ self.x ]
self.y_trajectory = [ self.y ]
[docs] def move(self, dt):
"""
Calculate the position of the ball after some has passed.
The new position is calculated by
.. math ::
x(t + \Delta t) = x(t) + v_x(t) \Delta t + \\frac{1}{2} \\frac{F_x(t)}{m} (\Delta t)^2
y(t + \Delta t) = x(t) + v_y(t) \Delta t + \\frac{1}{2} \\frac{F_y(t)}{m} (\Delta t)^2.
The new position is saved in the properties of the ball.
Args:
dt (float): timestep :math:`\Delta t`
"""
self.x += self.vx*dt + 0.5*self.fx/self.mass*dt*dt
self.y += self.vy*dt + 0.5*self.fy/self.mass*dt*dt
[docs] def accelerate(self, dt):
"""
Calculate the velocity of the ball after some time has passed.
The new velocity is calculated by
.. math ::
v_x(t + \Delta t) = v_x(t) + \\frac{F_x(t)}{m} \Delta t
v_y(t + \Delta t) = v_y(t) + \\frac{F_y(t)}{m} \Delta t.
The new velocity is saved in the properties of the ball.
Args:
dt (float): timestep :math:`\Delta t`
"""
self.vx += self.fx/self.mass*dt
self.vy += self.fy/self.mass*dt
[docs] def reset_forces(self):
"""
Forget the forces from a previous time.
"""
self.fx = 0.0
self.fy = 0.0
[docs] def apply_gravity(self, g):
"""
Make gravity act on the ball.
Adds the force :math:`F_y = - m g`.
Args:
g (float): gravitational acceleration :math:`g`
"""
self.fy += -self.mass * g
[docs] def apply_air_resistance(self, drag):
"""
Make air resistance act on the ball.
Adds the force :math:`\\vec{F} = - \\gamma m \\vec{v}`.
Args:
drag (float): air drag coefficient :math:`\\gamma`
"""
self.fx -= drag * self.mass * self.vx
self.fy -= drag * self.mass * self.vy
[docs] def apply_bounce(self, restitution):
"""
Make the ball bounce from the ground (ground level is at y = 0).
Args:
restitution (float): Coefficient for the elasticity of the bounce.
1 = completely elastic, 0 = completely inelastic.
"""
if self.y < self.radius:
self.vy = -restitution*self.vy
self.y = -self.y + 2*self.radius
[docs] def record_position(self):
"""
Save the position of the ball in memory.
"""
self.x_trajectory.append( self.x )
self.y_trajectory.append( self.y )
[docs]def run_simulation(ball, g, drag, restitution, dt, simulation_time, recording_dt):
"""
Run a dynamic simulation.
The simulation runs repeating the following steps:
* Forces are applied on the ball using e.g. :meth:`Ball.apply_gravity`.
* The ball is moved with :meth:`Ball.move`.
* The velocity is updated with :meth:`Ball.accelerate`.
The trajectory is saved in the :class:`Ball` object.
.. note ::
This function is incomplete!
Args:
ball (Ball): the moving ball
g (float): acceleration due to gravity
drag (float): air resistance factor (between 0 and 1)
restitution (float): impact elasticity factor (between 0 and 1)
dt (float): time step :math:`\\Delta t`
simulation_time (float): total simulation time
recording_dt (float): time between trajectory recording
"""
steps = int(simulation_time / dt)
recording_interval = int(recording_dt / dt)
for step in range(1,steps):
#
# ToDo: implement the simulation
#
ball.reset_forces()
# record the position every so often
if step%recording_interval == 0:
ball.record_position()
[docs]def main(v_start, angle):
"""
Main program.
Creates a ball, throws it at an angle and simulates the throw.
Args:
v_start (float): initial speed
angle (float): throwing angle in degrees, measured from the horizontal
"""
# starting position
x_start = 0.0
y_start = 1.5
# starting velocity
degree = pi/180
throw_angle = angle*degree
vx_start = v_start*cos(throw_angle)
vy_start = v_start*sin(throw_angle)
# physical constants:
# changing these will change the physics of the system
# gravity
g = 9.8
# air resistance (between 0 and 1)
drag = 0.2
# bounce elasticity (between 0 and 1)
restitution = 0.7
# simulation time step
dt = 0.01
# create the ball
ball = Ball(x_start, y_start, vx_start, vy_start)
# how often do you want to save data for plotting?
recording_dt = 0.1
# how long do you want to simulate?
simulation_time = 20.0
# run the simulation
run_simulation(ball, g, drag, restitution, dt, simulation_time, recording_dt)
# visualize the simulation
draw_trajectory(ball)
draw_motion_blur_image(ball)
animate(ball)
# Run this code if the file is executed as a program
# but not if it's loaded as a module.
if __name__ == "__main__":
v = 8.0
angle = 75 # degrees
main(v, angle)