function P1()
close all;

% Four wind modes:
wind1 = @(t)0;
wind2 = @(t)-10;
wind3 = @(t)10*(1-mod(floor(t),2));
wind4 = @(t)10*randn;

global fCount;
% size of figures:
xmax = 210;
ymax = 110;

% wind = 0:
figure(1);
hold;
max_downwind = 0;
for k = 1:17
    theta = k*pi/36;
    fCount = 0;
    [t trajectory] = cannonball_trajectory(theta, wind1);
    plot(trajectory(:,1),trajectory(:,2));
    if trajectory(end,1)>max_downwind
        max_downwind = trajectory(end,1);
        max_downwind_theta = theta;
        max_downwind_flight_time = t(end);
        max_downwind_velocity = trajectory(end, 4);
        max_fCount = fCount;
    end
end
axis([0 xmax 0 ymax]);
title(sprintf('No wind\n theta=%f, time=%f,\n distance=%f,\n velocity=%f, f-evals=%d',...
    max_downwind_theta, max_downwind_flight_time,...
    max_downwind, max_downwind_velocity, max_fCount));

% wind = -10 m/s:
figure(2);
hold;
max_downwind = 0;
for k = 1:17
    theta = k*pi/36;
    fCount = 0;
    [t trajectory] = cannonball_trajectory(theta, wind2);
    plot(trajectory(:,1),trajectory(:,2));
    if trajectory(end,1)>max_downwind
        max_downwind = trajectory(end,1);
        max_downwind_theta = theta;
        max_downwind_flight_time = t(end);
        max_downwind_velocity = trajectory(end, 4);
         max_fCount = fCount;
    end
end
axis([0 xmax 0 ymax]);
title(sprintf('Steady head wind\n theta=%f, time=%f,\n distance=%f,\n velocity=%f, f-evals=%d',...
    max_downwind_theta, max_downwind_flight_time,...
    max_downwind, max_downwind_velocity, max_fCount));

% wind = 10 m/s:
figure(3);
hold;
max_downwind = 0;
for k = 1:17
    theta = k*pi/36;
    fCount = 0;
    [t trajectory] = cannonball_trajectory(theta, wind3);
    plot(trajectory(:,1),trajectory(:,2));
    if trajectory(end,1)>max_downwind
        max_downwind = trajectory(end,1);
        max_downwind_theta = theta;
        max_downwind_flight_time = t(end);
        max_downwind_velocity = trajectory(end, 4);
        max_fCount = fCount;
    end
end
axis([0 xmax 0 ymax]);
title(sprintf('Intermittent tailwind\n theta=%f, time=%f,\n distance=%f,\n velocity=%f, f-evals=%d',...
    max_downwind_theta, max_downwind_flight_time,...
    max_downwind, max_downwind_velocity, max_fCount));

% wind = 10*randn m/s:
figure(4);
hold;
max_downwind = 0;
for k = 1:17
    theta = k*pi/36;
    fCount = 0;
    [t trajectory] = cannonball_trajectory(theta, wind4);
    plot(trajectory(:,1),trajectory(:,2));
    if trajectory(end,1)>max_downwind
        max_downwind = trajectory(end,1);
        max_downwind_theta = theta;
        max_downwind_flight_time = t(end);
        max_downwind_velocity = trajectory(end, 4);
        max_fCount = fCount;
    end
end
axis([0 xmax 0 ymax]);
title(sprintf('Gusty wind\n theta=%f, time=%f,\n distance=%f,\n velocity=%f, f-evals=%d',...
    max_downwind_theta, max_downwind_flight_time,...
    max_downwind, max_downwind_velocity, max_fCount));


function [t trajectory nsteps] = cannonball_trajectory(theta, wind)
% theta: initial angle
% wind: wind conditions, (wt = wind(t))
% t: time when ode is computed
% trajector: [x y theta v]
% nsteps: the number of steps in the ode
opts = odeset('events', @g);
v0 = 50;
y0 = [0; 0; theta; v0];
% tra
[t trajectory] = ode45(@f, [0 Inf], y0, opts, wind);
nsteps = length(t); %How to compute nfevals?

function  dydt = f(t, y, wind)
% t: time
% y: [x y theta v]
% dydt: dy/dt
dydt = zeros(4,1);
g = 9.81; % acceleration of gravity
m = 15; % mass of the cannonball
dydt(1) = y(4)*cos(y(3));
dydt(2) = y(4)*sin(y(3));
dydt(3) = -g/y(4)*cos(y(3));
dydt(4) = -drag(t, dydt(1), dydt(2), wind)/m - g*sin(y(3));
global fCount;
fCount = fCount+1;

function D = drag(t, xprime, yprime, wind)
c = 0.2;
rho = 1.29;
s = 0.25;
D = c*rho*s/2*((xprime-wind(t)).^2+yprime.^2);

function [gstop, isterminal, direction] = g(t, y, wind)
% t: time
% y: [x y theta v]
gstop = y(2);
isterminal = 1;
direction = -1;


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 2 points for an event function that works
% 1 point for the correct air time
% 1 point for correct landing velocity
% 1 point for correct f-evals
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%