In [1]:
# Configure Jupyter so figures appear in the notebook
%matplotlib inline
# Configure Jupyter to display the assigned value after an assignment
%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'
# import functions from the modsim.py module
from modsim import *
A Vector object represents a vector quantity. In the context of mechanics, vector quantities include position, velocity, acceleration, and force, all of which might be in 2D or 3D.
You can define a Vector object without units, but if it represents a physical quantity, you will often want to attach units to it.
I'll start by grabbing the units we'll need.
In [2]:
m = UNITS.meter
s = UNITS.second
kg = UNITS.kilogram
Here's a two dimensional Vector in meters.
In [3]:
A = Vector(3, 4) * m
We can access the elements by name.
In [4]:
A.x
In [5]:
A.y
The magnitude is the length of the vector.
In [6]:
A.mag
The angle is the number of radians between the vector and the positive x axis.
In [7]:
A.angle
If we make another Vector with the same units,
In [8]:
B = Vector(1, 2) * m
We can add Vector objects like this
In [9]:
A + B
And subtract like this:
In [10]:
A - B
We can compute the Euclidean distance between two Vectors.
In [11]:
A.dist(B)
And the difference in angle
In [12]:
A.diff_angle(B)
If we are given the magnitude and angle of a vector, what we have is the representation of the vector in polar coordinates.
In [13]:
mag = A.mag
angle = A.angle
We can use pol2cart to convert from polar to Cartesian coordinates, and then use the Cartesian coordinates to make a Vector object.
In this example, the Vector we get should have the same components as A.
In [14]:
x, y = pol2cart(angle, mag)
Vector(x, y)
Another way to represent the direction of A is a unit vector, which is a vector with magnitude 1 that points in the same direction as A. You can compute a unit vector by dividing a vector by its magnitude:
In [15]:
A / A.mag
Or by using the hat function, so named because unit vectors are conventionally decorated with a hat, like this: $\hat{A}$:
In [16]:
A.hat()
Exercise: Create a Vector named a_grav that represents acceleration due to gravity, with x component 0 and y component $-9.8$ meters / second$^2$.
In [17]:
# Solution goes here
In [18]:
degree = UNITS.degree
radian = UNITS.radian
If you have an angle in degrees,
In [19]:
angle = 45 * degree
angle
You can convert to radians.
In [20]:
angle_rad = angle.to(radian)
If it's already in radians, to does the right thing.
In [21]:
angle_rad.to(radian)
You can also convert from radians to degrees.
In [22]:
angle_rad.to(degree)
As an alterative, you can use np.deg2rad, which works with Pint quantities, but it also works with simple numbers and NumPy arrays:
In [23]:
np.deg2rad(angle)
Exercise: Create a Vector named a_force that represents acceleration due to a force of 0.5 Newton applied to an object with mass 0.3 kilograms, in a direction 45 degrees up from the positive x-axis.
Add a_force to a_grav from the previous exercise. If that addition succeeds, that means that the units are compatible. Confirm that the total acceleration seems to make sense.
In [24]:
# Solution goes here
In [25]:
# Solution goes here
Here's a Params object that contains parameters for the flight of a baseball.
In [26]:
t_end = 10 * s
dt = t_end / 100
params = Params(x = 0 * m,
y = 1 * m,
g = 9.8 * m/s**2,
mass = 145e-3 * kg,
diameter = 73e-3 * m,
rho = 1.2 * kg/m**3,
C_d = 0.33,
angle = 45 * degree,
velocity = 40 * m / s,
t_end=t_end, dt=dt)
And here's the function that uses the Params object to make a System object.
In [27]:
def make_system(params):
"""Make a system object.
params: Params object with angle, velocity, x, y,
diameter, duration, g, mass, rho, and C_d
returns: System object
"""
angle, velocity = params.angle, params.velocity
# convert angle to degrees
theta = np.deg2rad(angle)
# compute x and y components of velocity
vx, vy = pol2cart(theta, velocity)
# make the initial state
R = Vector(params.x, params.y)
V = Vector(vx, vy)
init = State(R=R, V=V)
# compute area from diameter
diameter = params.diameter
area = np.pi * (diameter/2)**2
return System(params, init=init, area=area)
Here's how we use it:
In [28]:
system = make_system(params)
Here's a function that computes drag force using vectors:
In [29]:
def drag_force(V, system):
"""Computes drag force in the opposite direction of `v`.
V: velocity Vector
system: System object with rho, C_d, area
returns: Vector drag force
"""
rho, C_d, area = system.rho, system.C_d, system.area
mag = rho * V.mag**2 * C_d * area / 2
direction = -V.hat()
f_drag = direction * mag
return f_drag
We can test it like this.
In [30]:
V_test = Vector(10, 10) * m/s
drag_force(V_test, system)
Here's the slope function that computes acceleration due to gravity and drag.
In [31]:
def slope_func(state, t, system):
"""Computes derivatives of the state variables.
state: State (x, y, x velocity, y velocity)
t: time
system: System object with g, rho, C_d, area, mass
returns: sequence (vx, vy, ax, ay)
"""
R, V = state
mass, g = system.mass, system.g
a_drag = drag_force(V, system) / mass
a_grav = Vector(0, -g)
A = a_grav + a_drag
return V, A
Always test the slope function with the initial conditions.
In [32]:
slope_func(system.init, 0, system)
We can use an event function to stop the simulation when the ball hits the ground:
In [33]:
def event_func(state, t, system):
"""Stop when the y coordinate is 0.
state: State object
t: time
system: System object
returns: y coordinate
"""
R, V = state
return R.y
In [34]:
event_func(system.init, 0, system)
Now we can call run_ode_solver
In [35]:
results, details = run_ode_solver(system, slope_func, events=event_func)
details
The final label tells us the flight time.
In [36]:
flight_time = get_last_label(results) * s
The final value of x tells us the how far the ball landed from home plate:
In [37]:
R_final = get_last_value(results.R)
x_dist = R_final.x
In [38]:
xs = results.R.extract('x')
ys = results.R.extract('y')
xs.plot()
ys.plot()
decorate(xlabel='Time (s)',
ylabel='Position (m)')
savefig('figs/chap22-fig01.pdf')
We can plot the velocities the same way.
In [39]:
vx = results.V.extract('x')
vy = results.V.extract('y')
vx.plot(label='vx')
vy.plot(label='vy')
decorate(xlabel='Time (s)',
ylabel='Velocity (m/s)')
The x velocity slows down due to drag.
The y velocity drops quickly while drag and gravity are in the same direction, then more slowly after the ball starts to fall.
Another way to visualize the results is to plot y versus x. The result is the trajectory of the ball through its plane of motion.
In [40]:
def plot_trajectory(results):
xs = results.R.extract('x')
ys = results.R.extract('y')
plot(xs, ys, color='C2', label='trajectory')
decorate(xlabel='x position (m)',
ylabel='y position (m)')
plot_trajectory(results)
savefig('figs/chap22-fig02.pdf')
The draw function should take as parameters a State object and the time. It should draw a single frame of the animation.
Inside the draw function, you almost always have to call set_xlim and set_ylim. Otherwise matplotlib auto-scales the axes, which is usually not what you want.
In [41]:
xs = results.R.extract('x')
ys = results.R.extract('y')
def draw_func(state, t):
set_xlim(xs)
set_ylim(ys)
x, y = state.R
plot(x, y, 'bo')
decorate(xlabel='x position (m)',
ylabel='y position (m)')
In [42]:
animate(results, draw_func)
Exercise: Delete the lines that set the x and y axes (or comment them out) and see what the animation does.
In [43]:
V = Vector(3, 4)
type(V)
A ModSimVector is a specialized kind of Pint Quantity.
In [44]:
isinstance(V, Quantity)
There's one gotcha you might run into with Vectors and Quantities. If you multiply a ModSimVector and a Quantity, you get a ModSimVector:
In [45]:
V1 = V * m
In [46]:
type(V1)
But if you multiply a Quantity and a Vector, you get a Quantity:
In [47]:
V2 = m * V
In [48]:
type(V2)
With a ModSimVector you can get the coordinates using dot notation, as well as mag, mag2, and angle:
In [49]:
V1.x, V1.y, V1.mag, V1.angle
With a Quantity, you can't. But you can use indexing to get the coordinates:
In [50]:
V2[0], V2[1]
And you can use vector functions to get the magnitude and angle.
In [51]:
vector_mag(V2), vector_angle(V2)
And often you can avoid the whole issue by doing the multiplication with the ModSimVector on the left.
Exercise: Run the simulation with and without air resistance. How wrong would we be if we ignored drag?
In [52]:
# Hint
system_no_drag = System(system, C_d=0)
In [53]:
# Solution goes here
In [54]:
# Solution goes here
In [55]:
# Solution goes here
In [56]:
# Solution goes here
In [57]:
# Solution goes here
Exercise: The baseball stadium in Denver, Colorado is 1,580 meters above sea level, where the density of air is about 1.0 kg / meter$^3$. How much farther would a ball hit with the same velocity and launch angle travel?
In [58]:
# Hint
system2 = System(system, rho=1.0*kg/m**3)
In [59]:
# Solution goes here
In [60]:
# Solution goes here
Exercise: The model so far is based on the assumption that coefficient of drag does not depend on velocity, but in reality it does. The following figure, from Adair, The Physics of Baseball, shows coefficient of drag as a function of velocity.
I used an online graph digitizer to extract the data and save it in a CSV file. Here's how we can read it:
In [61]:
baseball_drag = pd.read_csv('data/baseball_drag.csv')
mph = Quantity(baseball_drag['Velocity in mph'], UNITS.mph)
mps = mph.to(m/s)
baseball_drag.index = magnitude(mps)
baseball_drag.index.name = 'Velocity in meters per second'
baseball_drag
Modify the model to include the dependence of C_d on velocity, and see how much it affects the results. Hint: use interpolate.
In [62]:
# Solution goes here
In [63]:
# Solution goes here
In [64]:
C_d = drag_interp(43 * m / s)
In [65]:
# Solution goes here
In [66]:
# Solution goes here
In [67]:
# Solution goes here
In [68]:
# Solution goes here
In [69]:
# Solution goes here
In [70]:
# Solution goes here
In [71]:
# Solution goes here
In [72]:
# Solution goes here
In [ ]: