% ### EXprojectileMOD.m ###      2020.01.28  C. Bergevin

% [REF: ex.4.3.2 from Fowles & Cassidy 2005]
% Purpose: Modified version of EXprojectile.m to show "optimal" angle for
% projectile flights w/ drag is less than 45 degrees

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

clear
% ------------------------------------------
P.thetaA= linspace(36,42,21);  % launch angle [degrees] 
P.g= 9.8;     % grav. const. [m^2/s] {9.8}
P.drag= 1;     % 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 15];  % 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.y0(1)= P.coord0(1); P.y0(3)= P.coord0(2);  % initial horiz. and vert. positions    

% --- set up fig. plus color-coding scheme
colormap(jet(numel(P.thetaA)))
jetcustom = jet(numel(P.thetaA));
figure(1); clf; hold on; grid on;

% ---
for nn=1:numel(P.thetaA)
    P.theta= P.thetaA(nn);
    P.theta= pi*P.theta/180;   % convert launch angle to rads
    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);
    indxG= find(vals(:,3)<0,1);  % find when object hits the ground
    % --- rename vars. (excluding those in the ground!) & plot
    x= vals(1:indxG,1); xdot= vals(1:indxG,2);
    z= vals(1:indxG,3); zdot= vals(1:indxG,4);
    plot(x,z,'-','LineWidth',1,'Color',jetcustom(nn,:));   
end

xlabel('x [m]');  ylabel('z [m]');
hC=colorbar; caxis([min(P.thetaA) max(P.thetaA)]);
ylabel(hC, 'Launch angle [deg]')
title('Projectile range w/ quadratic drag for multiple launch angles')