% ### EXprojectile.m ###      2020.01.28  C. Bergevin

% [REF: ex.4.3.2 from Fowles & Cassidy 2005]
% Purpose: Solve/plot 2-D projectile motion for spherical object w/ (optional)
% quadratic drag (x is horiz. position, z vert. pos.) 

% ---- Notes
% o v0= 143.2 mph ~ 64 m/s

clear
% ------------------------------------------
P.g= 9.8;     % grav. const. [m^2/s] {9.8}
P.drag= 0;     % boolean to incl. drag: 0=no drag, 1=drag {1}
P.v0= 64;     % launch velocity [m/s] {64}
P.theta= 39;     % launch angle [degrees] {45}
P.D= 0.073;   % diameter of object [m] {0.073}
P.m= 0.145;   % mass of object [kg] {0.145}
P.coord0= [0 0];  % initial [x z] coords [m] {[0 0]}
P.tLim= [0 10];  % time limits of integration [s]
P.tRez= 300;     % # of (interp.) time points for integration interval {300?}
% ------------------------------------------

% --- derived params. 
if (P.drag==0), P.gamma= 0;   % determine assoc. const. from input params.
else P.gamma= 0.15*P.D^2/P.m;  end
P.theta= pi*P.theta/180;   % convert launch angle to rads
P.y0(1)= P.coord0(1); P.y0(3)= P.coord0(2);  % initial horiz. and vert. positions    
P.y0(2)= P.v0*cos(P.theta);  % initial horiz. velocity
P.y0(4)= P.v0*sin(P.theta);  % initial vert. velocity

% --- use built-in solver ode45 to numerically integrate
[t vals] = ode45('PROJECTILEfunction',linspace(P.tLim(1),P.tLim(2),P.tRez),P.y0,[],P);

% --- kludge: find when object hits the ground (and indicate if it hasn't)
indxG= find(vals(:,3)<0,1); flag= 0;
if (isempty(indxG)), disp('Longer int. time needed (to hit ground)'); indxG=size(vals,1); flag=1; end
indxH= find(vals(:,4)<0,1);  % index where velocity flips sign

% --- rename vars. (excluding those in the ground!)
x= vals(1:indxG,1); xdot= vals(1:indxG,2);
z= vals(1:indxG,3); zdot= vals(1:indxG,4);

% --- spit back a few vals. to screen
if (flag==0), disp(['total flight time= ',num2str(t(indxG)),' s']);
    disp(['horizontal dist. covered= ',num2str(x(indxG)),' m']);
    disp(['max. vert. height= ',num2str(z(indxH)),' m']);   end

% ---- visualize
figure(1); clf; h1= plot(x,z,'k-','LineWidth',1); hold on; grid on;
xlabel('x [m]');  ylabel('z [m]');