As usual, this is an open-ended project. Use your instincts and judgement to choose and ignore options with the goal of maximizing results while minimizing effort.
Do either:
The expectation is that you'll be using various sources for your investigations, so make sure you cite all references you used: the more (that you actually use) the better. Technique, results, and presentation (especially graphics) are all important. The elegance and organization (and correctness) of the code count, but not as much as the report, which should contain enough information to convince the reader that your code must be correct.
In all this work, you will use Matlab's ODE solvers: let's just stick to ode23. The ODE solution itself is a simple and well-defined problem: one 'derivative' function just encodes the problem's differential equations (often written out for you in mathematics). The serious computing is done by a 1-liner: calling ode23 with a few simple arguments. These ballistics problems call for various amounts of outside research (e.g. just how does atmospheric density vary with altitude?), and a little interesting programming (finding the range and range-maximizing elevation with no user interaction).
At least one good previous report on the Balls in the Air project is available on the 160 website .../u/brown/160/.
There are plenty of problems you can solve in closed form concerning the trajectory of bodies in a vacuum: Like These. But we normally fight in the atmosphere, and it turns out that air friction at muzzle-velocity speeds goes as the square of the speed of the the projectile (unlike the linear friction we saw in the mass-spring-damper example). It's also proportional to the cross-sectional area of the ball and the density of air. We get two or three coupled non-linear (because of the square of the velocity) second-order equations (for the vertical and horizontal, or also sideways, component to the ball's flight.)
State: Our first problem is to write the derivative functions for our problem variants. Since we want range and height of the cannonball (and later, its right-to-left deviation), let's say x, y, and z are the range, height, and left-deviation coordinates of the ball, respectively. Then x' (called u below), y' (called v), and z' (w) are the x, y, and z velocity components, respectively. The physics (coming up) will give us x'', y'', and z''. That means our state vector will consist of (x, x', y, y', z, z') (in some order). The derivative function takes in the state vector and returns its derivative (it copies x', y', z' into the output and uses the physics to compute x'', y'', z'').
Physics:
Let s be the ball's speed through the air; then
s2 =
(u2 + v2
+ w2).
A is the cross-sectional area of the ball and m its mass,
ρ is the density of air (zero for a vacuum), g the force of gravity,
and c is a dimensionless drag coefficient.
For spheres, c tends to be about .5.
Technically, c is the ratio between the actual drag force, and the
theoretical drag force that would result from a hypothetical object
that brought all the air it intercepted up to its own speed.
Unstreamlined objects like cubes, flat surfaces normal to the airflow, and
end-on cylinders tend to have drag coefficients pretty close to 1.0.
Nice, streamlined, torpedo-shaped objects can have drag coefficients
of 0.1 or even lower.
The force due to air resistance for a symmetrical object
is directly against the direction of motion
(that is, in the direction of the apparent "wind"), and has magnitude
given approximately by
F = (c A ρ s2) / 2.
Since we are formulating the problem in terms of x, y, and z components,
we have to determine the force components in those directions.
In particular, we want the force in each direction as a function of
the x, y, and z velocity components (x', y' and z').
In this section we assume the medium (air) is at rest, but below that assumption changes. For now let's stick with the "still air" assumption.
The unit vector in the direction of
instantaneous velocity is
(α, β, γ) = (1/√s2) (u,v,w).
Its three components give the proportion of the force to be allocated
to the (x,y,z) directions.
From " force equals mass times acceleration", or F = ma, we get a =
F/m,
thence
x'' = - (1/(2m)) (ρ c A s2 α) (drag alone)
y'' = -g - (1/(2m)) (ρ c A s2β) (gravity plus drag)
Express the area in terms of the diameter and pi:
A = π (d/2)2
Make all these substitutions and we have all we need to write the derivative function that computes x'' and y'' from x, x', y, y', and our physical constants. We'll deal with z later.
What Now? The cannon points up at some elevation angle θ. Let's say the cannon points out the x axis, y is up, z is left. Cannon has some muzzle (initial) velocity, which should be subsonic for these equations: let's standardize on 150 m/sec, which is just under 1/2 the speed of sound. Given the muzzle velocity and your chosen shooting elevation angle θ you can calculate the initial velocities x', y',z' at time zero (that is, u0, v0, w0). (hint: z' =0, and sin θ and cos θ are involved in the calculations for u0, v0). The initial [0 0 0] position and [u0 v0 0] velocity conditions are sent in as the two columns of the 3x2 initial-conditions matrix to ode23() (see below).
For the problems below unless otherwise directed, use these base specifications in SI units (used to be called MKS for meter-kg-sec): muzzle velocity = 150 m/sec, mass of cannonball = 5kg, ρ = 1.23 kg/m3, g = 9.81 m/sec/sec, c = .5. For an iron cannonball of mass 5kg, its diameter (we need that to compute A) follows from the density of iron: 7860 kg/m3.
We're going to use Matlab's ode23(), and it wants the handle of the derivative function, a time span, and an initial state vector (let's say the cannon is at x = y = z = 0). Some experimentation is needed on the timespan, since our equations don't know anything about the ground and we should stop our solution about when the ball hits the ground.
With constants close to those above and
an elevation of 45 degrees I get a trajectory that looks like this: plotting
height (in m) against range (in m).
Notice the slight deviation from a pure parabola.
In fact, you can just set air
density (ρ) to 0 and run that case with
>> hold on
to superimpose the plots to visualize the difference.
Soon you'll be asked to see how elevation,
drag, and mass affect the trajectory.
You'd expect a lighter, projectile to slow down faster
(from intuition and also from the equations),
so it would have a less symmetric trajectory.
.
Footnote:
Who's this kindly old gentleman?
(from wikimedia commons, Scanned from German "Meyer's Encyclopedia",
1906.)
Turns out for
modern bullets there's a special drag coefficient called the
The
Ballistic Coefficient . It is a number based on
the drag coefficient, only normalized to an "ideal bullet"
that this
benign duffer produced in 1881, with ballistic coefficient of 1 by
definition. (So now you know: if he's dictating standards, he's
not just somebody's beloved grandpa).
Unlike the drag coefficient, a smaller ballistic coefficient means more drag.
Regular bullets in use today have ballistic coefs varying from
about .01 (.177
pellets)
to about 1, but you can make a bullet with BC 1.1 out of special
materials with special lathes...better than the first "ideal bullet".
.
Battle of Crécy between the English and French in the Hundred Years' War. From a 15th-century illuminated manuscript of Jean Froissart's Chronicles (BNF, FR 2643, fol. 165v). (Commons.Wikimedia.org)
We're now in 3-D, and we are measuring the forces on the cannonball in
a medium that is moving at windvel, the wind velocity, here
80mph in the z direction.
We thus relativize the z velocity to the moving atmosphere, and
compute the squared velocity:
s2 = x'2 + y'2 + (windvel- z') 2.
The formula for z''
then looks like that for x'', only with a γ to project the force
onto the (relative) z' direction, instead of an α
(see the Physics section above for α, etc. defs.).
For debugging, Should try it with no x or y velocity and a very light cannonball (like a balloon).
FWIW, I get a pretty significant z deflection for the base specifications and 45 degree shooting elevation.
Optional: Code Organization, Approaches, Opinions, Hints:
function dy = cannon_der(t,y) script shoot.m function x0 = zero_cross(x,y) function x = range(angle, diameter, mass) function xmax = max_range(low_angle, high_angle) function dy = cannon_air_der(t,y)
function dy = cannon_der(t,y) This is the always-needed derivative equations whose handle you pass to your ODE solver. It specifies the derivatives for the state vector that consists, in some order that makes sense to you, of the quantities (x,y,z,x',y',z'). The output is a COLUMN 6-vector specifying their derivatives (x', y', z', x'', y'', z''). The z (East-West) direction is used in the last part of the question. Note that since there is no explicit time dependence in this problem, the t value will be ignored. It is needed only to fit the form used by ode23().
script shoot.m This is simply the call to ode23() and various plot commands to help get acquainted with the problem and for parts 1 and 5 of the assignment.
function x = range(der_func_handle, angle, diameter, mass) This is the basic cannon-shooting simulator that we can use in all parts of the assignment. It uses your ODE solver (ode23() is fine) to compute the range of the cannonball (where its height (y) goes to 0). Thus this function needs the handle to cannon_der, it needs to set initial conditions (using the angle) and a time span (experiment until you get one that gets y (height) into negative values for the range of angles you want to consider) and it calls zero-cross. Its parameters are the most commonly-changed ones (and you have to change them in the second part of the assignment).
IMPORTANT! cannon_der() needs to use these parameters (angle, diameter, mass) but we can't pass them into cannon_der() since its parameters have to look like (t,y). So declare these variables global in both cannon_der() and range().
function x0 = zero_cross(x,y) processes the output of ode23() to estimate the value of x for which the ball hits the ground (i.e. for which y(x) = 0). ode23's output consists of the varying state vector Y (each row is a [x y z x' y' z'] vector) and the associated x value vector, where Y = f(x). x0 is the computed estimate of the value of x for which y=0. One coarse estimate is the first x for which y changes sign from positive to negative.
Better is a linear interpolation: Let's say y starts out positive, as in this assignment. Given (x1,y1), with y1 the last positive y, and its neighboring point (x2,y2), with y2 the first negative y, you can draw a little picture with the two points, a line between them, and the line y=0, and notice you want to solve for x0, where a point on the line is (x0,0). A similar triangles argument and a little high school algebra gives you x0 in terms of x1, y1, x2, and y2.
function xmax = max_range(low_angle, high_angle) This answers the second part of the assignment. It calls range() for various angles (maybe let them vary freom low to high by 2 degrees?), gets the ranges and returns the maximum one. I'd use a single for-loop that initializes a variable max_range to the first range and best_angle to the first angle, and subsequently updates the max_range and best_angle to the current range and angle if the current range is longer than the previous max.
function dy = cannon_air_der(t,y) For the fourth part, we need a new derivative function that is just like cannon_der (globals and all) but instead of being a constant, the air-density is a function of ρ, the sea-level density, and y (the height). That change must be made in the three derivative functions for x', y', and z'. Come up with a close-enough function (e.g. google 'air density table'), justify it with a reference, and carry on!
The NYTimes has hired us to quantify (and thereby they hope to invalidate) their right-wing rival's story. Here's what we (well, you) will do. (Note: NYTimes needs references for all numbers, functions, facts, etc.). About 10 minutes googling should be all you need.
The strategy is to simulate the trajectories of hits, see if each hint is a homer, and to see if we can replicate the claimed statistical improvement in the number of homers given different velocities of following winds.
Hint: don't even think about writing your final program and starting to play with it. Small routines, thoroughly understood and debugged: that's the thing. For instance, you could test your monte carlo routine by only varying one parameter (elevation or velocity), using a simple uniform distribution, and making sure the data make sense). Then the other parameter. Only then do both. Continue this conservative approach.
All this is used in a Monte Carlo simulation to generate a realistic population of ball trajectories. Monte Carlo is just a fancy way of saying "repeated probabilistic simulation": A simpler example is the way π was used in A 160 Assignment (look down to "Π part 3: Monte Carlo Simulation").
You'll want to discard numbers below and above some thresholds (velocity should be between some reasonable minimum speed off the bat (40 mph?) and some known maximum whose value you can get from the web.
Do your results, given the limittions of your approximations, agree with the story?
If that link dies, here's a local copy. Offhand, what do you think of his science? Present your data on the downward angle and compare. How do the Web-based Expert's assumptions square with your findings?