File: ef230-2023-08/www/examples/projectile2003.m Download
%! Projectile motion with a target and drag - example from EF 102 2003
% University of Tennessee - Will Schleter

% Clear memory and command window.
clear all; clc; close all;
% Define conversion functions for degrees <--> radians
d2r = inline('a.*pi./180');
r2d = inline('a.*180./pi');
%======================================================================%
% Give choice of the target
ti = menu('Choose target','Simple building','Building with peak','Incline','Surprise');
switch ti
    case 1 % simple building
        xtarg = [30 30 50 50];
        ytarg = [ 0 30 30  0];
    case 2 % peak building
        xtarg = [30 30 40 50 50];
        ytarg = [ 0 30 40 30  0];
    case 3 % incline
        xtarg = [20 120 120];
        ytarg = [ 0  50  0];
    case 4 % surprise
        xtarg = [35 33 31 30 31 33 35 34 32 31 31 32 34 36 34 33 34 36 38 40 42 43 45 45 44 44 43 44 46 48 49 50 49 48 47 48 49 49 47 45];
        ytarg = [0 2 5 8 10 12 13 15 17 19 20 23 24 26 28 31 34 36 37 37 37 36 34 31 29 27 25 24 24 22 20 18 13 11 10 9 7 4 1 0];
    otherwise
        error('You did not select a target');
        return
end

% Start the plot and show the target
figure(1);
hold on;
axis equal;
% Plot the outline of the target.
patch(xtarg,ytarg,'g');

% Annotate plot.
title('Projectile Motion ');
xlabel('Horizontal Position - ft');
ylabel('Vertical Position - ft');

% initial values
v0=45;
ang=60;
x0=0;
y0=0;
drag=0;
dt=.01;
prompt={'Initial velocity (ft/sec)',
    'Angle (degrees)',
    'Initial X (ft)',
    'Initial Y (ft)',
    'Drag (1/sec)',
    'Time step (sec)' };

htxt=0; % init handle for text note displaying results

% loop forever or until user hits cancel in the input window
while 1
%======================================================================%
% Get user input for all values
defs = {num2str(v0),num2str(ang),num2str(x0),num2str(y0),num2str(drag),num2str(dt)};
answers = inputdlg(prompt,'Projectile Motion',1,defs);
if isempty(answers); return; end  % done, goodbye

v0 = eval(char(answers(1)));
ang = eval(char(answers(2)));
x0 = eval(char(answers(3)));
y0 = eval(char(answers(4)));
drag = eval(char(answers(5)));
dt = eval(char(answers(6)));

g = -32.2;

% Calculate initial velocity and acceleration components.
v_x0 = v0 .* cos(d2r(ang));
v_y0 = v0 .* sin(d2r(ang));
a_x = 0;
a_y = g;


% Initialize the starting time(t) and set the x and y coordinates at the
% start (i.e., at t = 0). Defining the starting values this way makes sure
% data from a previous run, if it exists, is discarded
ii=1;
t = 0.0;
xi = x0;
yi = y0;
v_xi = v_x0;
v_yi = v_y0;

% From here we need to calculate the x and y coordinates at "small"
% intervals along the curve. The coordinates will be stored in the
% vectors xi & yi. We will build xi & yi one element at a time.
while ~inpolygon(xi(ii),yi(ii),xtarg,ytarg) & yi(ii)>=0
    t(ii+1) = t(ii) + dt;
    v_xi(ii+1) = v_xi(ii) + (a_x - v_xi(ii).*drag).*dt ;
    v_yi(ii+1) = v_yi(ii) + (a_y - v_yi(ii).*drag).*dt ;
    xi(ii+1) = xi(ii) + (v_xi(ii) + v_xi(ii+1)).*dt./2;
    yi(ii+1) = yi(ii) + (v_yi(ii) + v_yi(ii+1)).*dt./2;
    ii=ii+1;
end
% Plot the projectile's curve.
% comet(xi,yi);  % animated
plot(xi,yi,'b-'); % permanent
% Plot a large red dot at the point of impact.
plot(xi(ii),yi(ii),'ro','MarkerSize',10,'MarkerFaceColor','r');

% Add text note to display initial conditions and final answers.
% Set default values for all text on this figure
set(gcf,'DefaulttextHorizontalAlignment','left','DefaulttextVerticalAlignment','top','DefaulttextUnits','normalized');

if htxt, delete(htxt), end; % delete the last note, if present
htxt = text(.1,.95,{['Launch height = ',num2str(y0,3),' ft'],...
        sprintf('Launch speed = %.1f ft/sec',v0),...
        sprintf('Launch angle = %.1f ^o',ang),...
        sprintf('Impact point = [%.1f %.1f] ft',xi(ii),yi(ii)),...
        sprintf('Impact time = %.2f sec',t(ii)),...
        sprintf('Time step = %.3f sec',dt),...
        sprintf('Drag = %.3f sec^-^1',drag),...
        sprintf('Impact velocity = %.1f ft/sec',norm([v_xi(ii) v_yi(ii)])),...
        sprintf('Net distance = %.1f ft',norm([xi(ii)-x0 yi(ii)-y0])) });

end % of while 1 loop
% End of program.