function distance = baseball(velocity, angle)
    theta = angle * pi / 180;
    vx = velocity * cos(theta);
    vy = velocity * sin(theta);

    options = odeset('Events', @events);
    [T, M] = ode45(@baseproj, [0, 100], [0, 1, vx, vy], options);
    plot(M(:,1), M(:,2))
    distance = M(end,1);
end

function [value,isterminal,direction] = events(t,W)
    % Locate the time when height passes through zero in a 
    % decreasing direction and stop integration.
    value = W(2);        % 
    isterminal = 1;      % Stop the integration
    direction = -1;      % Negative direction only
end

function dWdt = baseproj(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 = drag_force(V) / m;
    dVdt = gravity + drag;
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

