function  canonballbig
% Canon ball flight. Zero drag, constant gravity ballistics.
% Like the homework due Sept 14, 2006.
% OK, this one is easy with pencil and paper.
% Its just computer practice.
% Fancy stuff here: we end the integration when something interesting happens,
% namely when the projectile hits the ground.
% To do this we need to
%  1) use odeset to set options to include events
%             and also name the file where the 'event' is defined.
%     In this case the 'event' is going to be y passing through zero
%     from above. Canon ball hits ground. Stop calculating.
%  2) use a fancier calling of ode45 that returns not
%     just the value of t and z but also what happenned at the 'event'.
%       tevent = the time the event occured
%       zevent = the value of z (the state) when the event occured
%       ievent = which event occurred. In this example we only have
%                one event, so this is dumb.
%  3) Define an 'event' function.  Here we call it fstopevent.
%     has to return 3 things.
%            value = the thing which, when it passes through
%                    zero, should make the integration stop.
%                    In this case its the height y.
%       isterminal = set to 1 if you want the integration to stop, 
%                    and we do in this example.
%        direction = set to -1 if you only want to end the integration
%                    if y is decreasing when it passes through zero.
%                    Setting to zero will stop the integration at
%                    any time y=0, which could be right at the start.

% In another file is the same program, functionally, just without
% all the comments and stuff done to make it all more useable and
% readable.

%Set parameters used in the ODE solution
% Another fancy thing: we are going to pass the parameters
% to the ODE instead of defining  them inside the ODE function.
%  This adds parameters to the call of ode45 and to the definition
%   of the right hand side. These are the parameters
g = 10; c=2;

%Set initial conditions
v0= 100;                  % initial speed
theta = 45*pi/180;       % launch angle
vx0 = v0*cos(theta);
vy0 = v0*sin(theta);
x0=0; 
y0=0;
zzero=[x0 vx0 y0 vy0];   % pack initial conditions into z0

% Note the huge integration time below. Why?
% Because we are using 'events' to stop the integration when
% the canan ball hits the ground.  Thus  we need to leave enough
% time to make sure that happens. By using only 2 numbers for
% tspan we are letting the computer decide what times to give
% us output.
tspan=[0 1000];  %specify start and stop integration time

% Setup options for solving the ODE
% The new thing here is 'events'.  In the command that follows
% we are saying that the 'even' is in a function called fstopevent
% that is further down in this same text file.
options=odeset('abstol',1e-12,'reltol',1e-12,'events',@fstopevent);

% This next line is the core of the whole file.
% Note that the parameters g and c are passed to
% the right-hand-side function f.
%   SOLVE THE ODE
[t,z,tevent, zevent, ievent] =  ode45(@f,tspan,zzero,options,g,c);
% In this case, with a terminating 'event', z(end) = zevent.

%compare solution with the analytic solution
format long  % show lots of digits in printout
analyticsolution = 2*v0^2*sin(theta)*cos(theta)/g;% anaylic range
computersolution = zevent(1);  % also = z(end,1)
numericalerror= computersolution - analyticsolution  % This gets printed, no semicolon

%Unpack the solution so you can read it easily
x=z(:,1);
y=z(:,3);
plot(x,y)
title('Andy''s Canon ball flight 9/14/06: no air friction')
xlabel('x');  ylabel('y')
axis('equal')  % make x and y axis have same scale so shape looks right
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function defines the right hand side of the ODE
function zdot = f(t,z,g,c);
% Note, after the t and z is the same list of parameters
%  used in the calling ode 45: g and c.
% We don't use c in this example, but you might want
% to for your homework.

%Unpack z into variables one can understand
x  = z(1);
vx = z(2);
y  = z(3);
vy = z(4);

% calculate derivatives, this is the differenctial equation
xdot  = vx;
vxdot = 0;
ydot  = vy;
vydot = -g;

% pack up the derivatives into zdot.
zdot = [xdot vxdot ydot vydot]';
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% This function defines the 'events'
% In this case the 'event' is to stop integrating
% when the projectile hits the ground.
function [value,isterminal,direction] = fstopevent(t,z,g,c)
value = z(3) ; % look for zero crossings of y=z(3)
isterminal = 1; % 1 = stop, 0 = evaluate and move on
direction = -1; % -1 means decreasing, 0 = who cares, 1 = increasing
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%