Introductory Programming
Fall 2006

Baseball
--------

Here is a version of proj.m that includes drag:

function dWdt = proj(t, W)
    R = W(1:2);
    V = W(3:4);

    dRdt = V;                          
    dVdt = acceleration(t, R, V);

    dWdt = [dRdt; dVdt];
end

function dVdt = acceleration(t, R, V)
    g = 9.8;                       % acceleration of gravity in m/s^2
    gravity = [0; -g];

    m = 0.145;                     % mass in kilograms
    drag_acc = drag_force(V) / m;
    dVdt = gravity + drag_acc;
end

function Fd = drag_force(V)
    C = 0.3;      % dimensionless
    rho = 1.3;    % kg / m^3
    A = 0.0042;   % m^2
    v = norm(V);  % m/s
    
    Fd = - 1/2 * C * rho * A * v * V;
end

Notice that the entire computation is in vector form -- I never break
position or velocity up into x and y components.

The norm function computes the magnitude of a vector.

In notes15 I gave you a version of baseball that uses Events to
terminate ode45 when the ball hits the ground.  The alternative
looks like this:

function distance = baseball1(velocity, angle)
    % compute the initial conditions
    theta = angle * pi / 180;
    vx = velocity * cos(theta);
    vy = velocity * sin(theta);

    % call ode45 to compute the trajectory
    [T, M] = ode45(@proj, [0, 100], [0, 1, vx, vy]);
    
    % rifle through the results to see where the ball landed
    X = M(:,1);
    Y = M(:,2);
    for i=1:length(Y)
        if Y(i)<0
            distance = X(i);
            break
        end
    end
    plot(X, Y)
end

The function call structure for this program is

baseball1 -> ode45 -> proj -> acceleration -> drag_force



Evaluation 8
------------

After you have done Evaluation 8, try out these variations:

0) To see an animation of the trajectory, use comet3 instead of plot3

1) What happens if you add gravity (in the negative z direction)?

2) What happens if the strength of the magnetic field increases
   over time, so that the magnitude of B = t?

3) What happens if the strength of the magnetic field increases
   with distance from the origin, so the magnitude of B = norm(P)?

4) What happens if we add air resistance with a drag constant of 1,
   so drag_force = -norm(V) * V ?