% Longjump problem
clear all
clc
close all
warning off
disp('Use the following system of ODEs as model for long jump:')
disp('x''(t) = v(t)*cos theta(t)')
disp('y''(t) = v(t)*sin theta(t)')
disp('theta''(t) = -g/v(t) * cos theta(t)')
disp('v''(t) = -D/m -g*sin theta(t)')
disp('D = c*rho*s/2*(x''(t)^2 + y''(t)^2)')
disp('Common parameters:')
disp('g = 9.81 m/s^2 (acceleration gravity)')
disp('m = 80 kg (mass)')
disp('c = 0.72 (drag coefficient)')
disp('s = 0.50 m^2 (jumper''s cross-sectional area)')
disp('x(0) = y(0) = 0 (initial position)')
disp('theta(0) = pi/8 (takeoff angle)')
% Parameters for all parts
g = 9.81;   % m/s^2 acceleration gravity
m = 80;     % kg mass
c = 0.72;   % drag coefficient
s = 0.50;   % m^2 jumper's cross-sectional area
t0=0;
tmax=1;
% Use event handling to find when y-coordinate is zero
opts = odeset('events',@LongjumpEvent);
% (a) "Nominal" high altitude jump
disp(' ')
disp('(a) "Nominal" high altitude jump')
disp('rho = 0.94 kg/m^3 (air density)')
disp('v(0) = 10 m/s (initial velocity)')
rho = 0.94;
v0 = 10;
y0=[0 0 pi/8 v0];    % initial condition
[t,y,te,ye] = ode45(@Jump,[t0 tmax],y0,opts,g,m,c,s,rho);
figure('Position', [570 300 450 400])
hold on
plot(t,y(:,1),'-o');
plot(t,y(:,2),'-x');
title(sprintf('(a) Distance and height off ground for long jump problem')) 
hold off
pause
figure('Position', [570 300 450 400])
hold on
plot(t,y(:,3),'-o');
plot(t,y(:,4),'-x');
title(sprintf('(a) Angle and velocity')) 
hold off
% Find t when y is zero using event handling
tf = te(2);
xf = ye(2,1);
% Alternative solution
%pp = spline(t,y(:,2));
%tf = fzero(@ppval,tmax,[],pp);
%xf = ppval(spline(t,y(:,1)),tf);
disp(sprintf('Distance for nominal jump at high altitude: %3.2f m',xf))
pause
figure('Position', [570 300 450 400])
xt = ppval(spline(t,y(:,1)),linspace(t0,tf,100));
yt = ppval(spline(t,y(:,2)),linspace(t0,tf,100));
plot(xt,yt);
axis equal
title(sprintf('(a) Trajectory of the jump')) 
pause
% (b) "Nominal" sea level jump
disp(' ')
disp('(b) "Nominal" sea level jump')
disp('rho = 1.29 kg/m^3 (air density)')
disp('v(0) = 10 m/s (initial velocity)')
rho = 1.29;
y0=[0 0 pi/8 v0];    % initial condition
[t,y,te,ye] = ode45(@Jump,[t0 tmax],y0,opts,g,m,c,s,rho);
figure('Position', [570 300 450 400])
hold on
plot(t,y(:,1),'-o');
plot(t,y(:,2),'-x');
title(sprintf('(b) Distance and height off ground for long jump problem')) 
hold off
xf = ye(2,1);
% Alternative solution
%pp = spline(t,y(:,2));
%tf = fzero(@ppval,tmax,[],pp);
%xf = ppval(spline(t,y(:,1)),tf);
disp(sprintf('Distance for nominal jump at sea level: %3.2f m',xf))
pause
% (c) Sprinter's high altitude jump
disp(' ')
disp('(c) Sprinter''s high altitude jump')
disp('rho = 0.94 kg/m^3 (air density)')
disp('Determine v(0) so that distance equals Bob Beamon''s 8.90 m jump')
rho = 0.94;
v0 = 11;    % initial guess
while (xf < 8.9)
    v0 = v0+.005;
    y0=[0 0 pi/8 v0];    % initial condition
    %[t,y] = ode45(@Jump,[t0 tmax],y0,[],g,m,c,s,rho);
    [t,y,te,ye] = ode45(@Jump,[t0 tmax],y0,opts,g,m,c,s,rho);
    xf = ye(2,1);
    % Alternative solution
    %pp = spline(t,y(:,2));
    %tf = fzero(@ppval,tmax,[],pp);
    %xf = ppval(spline(t,y(:,1)),tf);
end
disp(sprintf('Initial velocity for Beamon''s jump: %f m/s',v0))
disp(sprintf('(computed distance: %f m)', xf))
figure('Position', [570 300 450 400])
hold on
plot(t,y(:,1),'-o');
plot(t,y(:,2),'-x');
title(sprintf('(c) Distance and height off ground for long jump problem')) 
hold off
pause
% (d) Sprinter's sea level jump
disp(' ')
disp('(d) Sprinter''s sea level jump')
disp('rho = 1.29 kg/m^3 (air density)')
disp(sprintf('v(0) = %f m/s (from previous problem)', v0))
rho = 1.29;
y0=[0 0 pi/8 v0];    % initial condition
[t,y,te,ye] = ode45(@Jump,[t0 tmax],y0,opts,g,m,c,s,rho);
figure('Position', [570 300 450 400])
hold on
plot(t,y(:,1),'-o');
plot(t,y(:,2),'-x');
title(sprintf('(d) Distance and height off ground for long jump problem')) 
hold off
xf = ye(2,1);
% Alternative solution
%pp = spline(t,y(:,2));
%tf = fzero(@ppval,tmax,[],pp);
%xf = ppval(spline(t,y(:,1)),tf);
disp(sprintf('Distance for jump at sea level with Beamon''s initial velocity: %3.2f m',xf))
warning on
