Baseball Simulation

% Find optimal initial velocity and angle at which to hit a baseball to
% make it over the Green Monster Wall
function baseball()
    % Loop through initial speed and angle combinations
    speed = 30 : 60;
    angle = 5 : 89;
    heights = zeros(length(speed), length(angle));
    for s = 1 : length(speed)
        for a = 1 : length(angle)
            % baseball simulation returns the final height of the ball
            h = baseball_simulation(speed(s), angle(a));
            heights(s, a) = h;
        end
    end
    figure(3)
    % plot the height of the ball as a function of angle and speed
    pcolor(angle, speed, heights);
    hold on
    % Plot line to show where the ball ends up close to the 12 m height is
    %(the wall) and 1 m, close to the ground
    clabel(contour(angle, speed, heights, [1, 12], 'LineColor', 'w', 'LineWidth', 1))
end

function res = baseball_simulation(init_speed, init_angle)
% Function takes in initial speed and angle of a baseball and returns the
% final height of the ball
    % Baseball constants
    b_mass = 0.145; % kg
    b_radius = 0.07468 / 2; % m
    b_c_area = pi * b_radius^2; % cross sectional area, m
    coeff_drag = 0.3;

    % Environment constants
    air_density = 1.2; % kg / m^3
    g = 9.8; % m / s^2
    wall_x = 97; % distance from place of hit
    wall_y = 12; % height of wall

    % Run the main function which should return the final height
    res = main();

    function res = main()
        res = simulate();
    end

    % Returns the final height of the ball
    function res = simulate()
        options = odeset('Events', @event_function);
        % Return all data points calculated by ode45, as well as the param
        % values at the time TE when the function event is called
        [Time, Y, TE, YE] = ode45(@change, [0, 10], ...
            [0, 1, velocity_vector(init_speed, init_angle)], options);
        res = YE(2); % final height of the ball, determined when the function stops

        % Trigger event when the ball hits the ground or reaches
        % the wall
        function [value, isterminal, direction] = event_function(t, params)
           % To trigger, the value should equal 0. If the x pos of ball is
           % the x pos of wall, it will be 0
           value = params(2);
           if wall_x - params(1) <= 0 || params(2) <= 0
               value = 0;
           % If the ball does not reach wall, check if it has hit the
           % ground
           end
           isterminal = 1; % stop function as soon as this event is reached
           direction = 0; % At any zero (when value = 0) consider this event, not
                          % only when the function is either increasing (1) or
                          % decreasing (-1)
        end
    end

    function res = velocity_vector(speed, angle)
        % Takes the magnitude of the velocity and the angle at which the
        % object is moving to the horizontal
        % Returns a velocity vector in the form (vx, vy)
        % Angle is in degrees, therefore must be converted to radians
        angle_rad = angle * pi / 180;
        res = [speed * cos(angle_rad), speed * sin(angle_rad)];
    end

    function res = change(t, params)
       % Position is a vector. (x, y)
       Pos = params(1 : 2);
       % velocity is a vector. (vx, vy)
       V = params(3 : 4);
       % change in position is velocity
       dPosdt = V;
       % change in velocity is acceleration
       dVdt = acceleration(t, Pos, V);
       res = [dPosdt; dVdt];
    end

    function res = acceleration(t, Pos, V)
        % Drag is dependent only on velocity, therefore position and time
        % are not needed. But just in case we included an acceleration
        % source dependent upon them, t and pos are passed into the
        % function.
       vx = V(1);
       vy = V(2);
       % Force of drag
       v = sqrt(vx^2 + vy^2);
       drag = -1/2 * coeff_drag * b_c_area * air_density * v^2;
       % Force of gravity
       gravity = -b_mass * g;

       % Unit velocity vector that points in the horizontal direction
       unit_vx = vx / v;
       % Unit velocity vector that points in the vertical direction
       unit_vy = vy / v;

       % Total horizontal acceleration is just the horizontal component of
       % drag. Must be divided by m because F = ma, a = F/m
       dvxdt = (drag * unit_vx) / b_mass;
       % Total vertical acceleration is the vertical component of drag and
       % the acceleration due to gravity
       dvydt = (gravity + drag * unit_vy) / b_mass;

       res = [dvxdt; dvydt];
    end

end