% ### 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]');