function tDevour = timeToDeath( angle )
% t_devour = timeToDeath( angle )
%       angle should be in degrees
%       Solves part (a) for problem 1, HW3, ACM11
% to find the best (or worse) angle to run at, you can do the following:
% [ang,val,flg] = fminbnd(@timeToDeath,-180,180)
%   (which is the fastest way to die), or
% [ang,val,flg] = fminbnd(@(t)-timeToDeath(t),-180,180)
%   (the best way to run)

% == The human starts at (0,0) and the raptors start at... ==
a = 10*sqrt(3)/3;
R1 = [-10,-a];
R2 = [10,-a];
R3 = [0,2*a];

% == Setup the human's path ==
theta = deg2rad(angle); c = [cos(theta),sin(theta)];
c = c/norm(c);
H_SPEED = 6;
h = @(t) H_SPEED*t*c;   % person runs is a straight line

% == Setup the raptors' equations of motion ==
R_SPEED = 25;
R_SPEED2 = 20;      % This is for the injured raptor
odefun =  @(t,y) R_SPEED * (h(t)-[y(1),y(2)])' ./ norm(h(t)-[y(1),y(2)]);
odefun2 = @(t,y) R_SPEED2* (h(t)-[y(1),y(2)])' ./ norm(h(t)-[y(1),y(2)]);
tol = .1;           % The human is caught when the raptors are this close
% When the raptor has caught the human, we'll stop modelling the true
% motion, because dividing by norm( x - y ) will result in huge numbers
% An alternative is to use the 'Events' option -- see "odeset".
cutoff = @(t,y) norm( h(t)-[y(1),y(2)])  > tol ;
odefun =  @(t,y) odefun(t,y) .* cutoff(t,y);
odefun2 = @(t,y) odefun2(t,y) .* cutoff(t,y);

% == Run the simulation ==
opts = odeset('abstol',1e-3);   % set the accuracy of the solver
% We need to tell ode45 the start and end time.  The start time can be 0,
% wlog.  The end time can be anything greater than t_devour -- but we don't
% know what d_devour is yet!  So, pick something conservative, like 10
% seconds.  Once you run this a few times, you can pick a more accurate
% estimate, like 1 second, or even .5 seconds.  Here, we use 1 second, so
% our time vector is [0,1].
[T1,Y1] = ode45(odefun ,[0,1],R1,opts); % raptor #1
[T2,Y2] = ode45(odefun ,[0,1],R2,opts); % raptor #2
[T3,Y3] = ode45(odefun2,[0,1],R3,opts); % raptor #3 (injured)

% == Calculate t_devour using the results ==
n = sqrt(sum( (h(T1)-Y1).^2 ,2));
f = find( n < tol );
t1 = T1(f(1));

n = sqrt(sum( (h(T2)-Y2).^2 ,2));
f = find( n < tol );
t2 = T2(f(1));

n = sqrt(sum( (h(T3)-Y3).^2 ,2));
f = find( n < tol );
t3 = T3(f(1));

tDevour = min( [t1,t2,t3] );

if nargout < 1
    fprintf('Time until raptors eat you: %.4f seconds\n',tDevour );
end

Published with MATLAB® 7.4