function [T,sol]=raptorchase(theta)
% RAPTORCHASE Simulate the survival time of the human in the XKCD raptor problem
% 
% The XKCD raptor problem is posed as follows:
%
%   A human is at the center of an equilaterial triangle with side length
%   20 meters.  At each corner is a velociraptor, whose maximum speed 
%   is 25 m/s.  The top corner has a velociraptor with a broken leg,
%   which is limited to 10 m/s speed.  What direction should you run in
%   order to maximize survival time?
%
% This function uses ode45 to explicitly simulate the raptor motion and
% compute the first time when a raptor is within 20 cm of the position of
% the human.  
%
% T = raptorchase(theta) returns the survival time for the angle theta.
% [T,sol] = raptorchase(theta) also returns the solution structure from
% ODE45 to enable plotting the answers.
%
% This function is based on ideas presented by Nick Henderson at the ICME
% Open Day in May.

vhuman = 6; % human velocity in m/s
vraptor0 = 10; % slow raptor velociy in m/s
vraptor = 15; % raptor velocity in m/s

triangle_dimension = 20; % in meters

raptor_min_distance = 0.2; % a raptor within 20 cm can attack

tmax = 10; % maximum time for integration

% Add the ODE function
    function dpos = change_in_positions(~,pos)
        human = [pos(1) pos(2)];
        nraptors = length(pos)/2 - 1;
        assert(ceil(nraptors) == nraptors) % integer operations are exact
        dpos = zeros(size(pos));
        dpos(1) = vhuman*cos(theta);
        dpos(2) = vhuman*sin(theta);
               
        for i=1:nraptors
            if i>1, vrap=vraptor;
            else    vrap=vraptor0;
            end
            raptor = [pos(2*i+1) pos(2*i+2)];
            tdir = (human - raptor)/norm(human-raptor);
            dpos(2*i+1) = tdir(1)*vrap;
            dpos(2*i+2) = tdir(2)*vrap;
        end
    end

% Add a function to stop the ODE evaluation when the raptors get close
    function [val,isterm,dir]=eaten_event(~,pos)
        dir=0; % matlab bookkeeping
        isterm=1;
        
        human = [pos(1) pos(2)];
        nraptors = length(pos)/2 - 1;
        assert(ceil(nraptors) == nraptors) % integer operations are exact
        raptor_dist = zeros(nraptors,1);
               
        for i=1:nraptors
            raptor = [pos(2*i+1) pos(2*i+2)];
            raptor_dist(i) = norm(human-raptor);
        end
        val = min(raptor_dist) - raptor_min_distance; % val = 0 when captured
    end

pos = 20*[zeros(2,1) simplex_coordinates1(2)];
p0 = pos(:); % unroll pos   
    
opts = odeset('Events',@eaten_event);
sol = ode45(@change_in_positions, [0,tmax], p0, opts);
T = max(sol.x);
end