ME 240 Winter 1998
Introduction to Dynamics

Due dates: Feb. 5, 1998 for Prof. Karnopp's section, Feb. 6, 1998 for Prof. Grosh's section

Computer Assignment #1: Projectile Flight

Problem Statement: Using Matlab plotting packages, plot the x-y trajectory of a projectile (mass) flying through the air. Assume that only a gravity force acts on the mass; and that initial velocity and particle direction (attack angle) are given. By plotting trajectories for multiple angles and initial velocities determine the following:

1. Using 30 m/s initial velocity, determine the angle which yields the longest distance in the horizontal plane (x-direction distance). On a single page, print out the x-y trajectory for various angles.
2. Does this best angle depend on initial velocity? On another figure, plot the x-y trajectory for one other initial velocity.
3. Now use the Matlab plot commands to plot the x and y dependence on time (i.e., plot x(t) and y(t)). Hand in a plot for a 30 degree initial trajectory.
4. Print out an "m-file" which performs the tasks above.
5. To determine the answers to these questions, did we really need to use a computer simulation or could we have obtained the answers another way (if so how!)?
6. Using the working model demo, ``experimentally'' determine the angle for which the ball travels the farthest in the x direction just before striking the ground for a fixed initial speed of the ball ( use 8.00 m/s ).  Does this agree with your previous results? How do mass and initial velocity affect the maximum horizontal distance?

Step by Step Matlab Help:

1. Running Matlab: The process for starting Matlab depends on which platform you are using. Click here for the Matlab Insructions protocols for the three CAEN platforms. Keep the web window up while you are running Matlab so that you can step through parts of this exercise. If you have never run Matlab before, Click here for a very short example (highly recommended).
2. First set up a few preliminaries constants. Note that everything following a % sign is ignored by Matlab (therefore it is just a comment on each line).

g = 9.81;                    %   gravitational acceleration (m/s^2)
N = 1000;                    %   No. of data points; sufficient for this problem
t_max = 10.0;                %   final time of the flight; OK for the 30 m/s v_0.
v0 =     30.0                %   Initial velocity of the ball (m/s) --> you can change
theta = YOU CHOOSE!!         %   Angle of attack = pi*(degrees)/180, i.e. in radians

v0 and theta determine the magnitude and direction of the initial velocity vector. Below, the method by which Matlab can compute the x and y dependence versus time is given. You will need to determine (find in the book!) the x and y formulas.
3. Now set up the arrays necessary for computation of the variables:

%
%
%   Calculations of the variables
%


t = 0:t_max/(N-1):t_max;      %  builds an array for the time of the flight 

x = YOU NEED TO PROVIDE THIS FORMULA; %  horizontal displacement
y = YOU NEED TO PROVIDE THIS FORMULA; %  vertical displacement


In these formulas, the y dependence will have a quadratic dependence on the time.
Matlab allows you to  do this in one step (avoiding a loop) if
instead of just using

t^2  % DON't USE THIS
 
you use the command

t.^2    %USE THIS

 which creates a new vector of the square of the time at each time step.
 So if you wanted to create a vector z(t) = a*t + b*t^2 (as you WILL want
 to do).  The Matlab command is

 z = a*t + b*t.^2
 
 No loops!

4. Now that the x and y dependence is know, we can plot the variables.

plot(x,y,'o');

5. The only problem here is that it is hard to find the exact zero crossing. So we can use a special Matlab command (any necessary tricks will be given beforehand). The values of x and y are stored in an array, in this case there are 1000 points in each of these arrays. To access points 100 to 200 of the array you use the command x(100:200) for instance. We want to know when the y value is zero (when the mass strikes the ground). The following command creates and array of the indexes for which the ball is below the ground

iz = find(y<0); % creates an array of all points in y <0

We only want the first such point

iz1 = iz(1);


and now we can make a better plot

plot(x(1:iz1),y(1:iz1),'o');

NOTE: the number of points in the x and y plotted arrays must be 
the same.  You can find an approximate time when the ball hits
the ground by printing the value of t(iz1).

%
%   Display the results :
%   max. horizontal distance & max. vertical height
% 
%

6. Now you can plot the results. The following prints a few pertanent items to the screen and creates a figure with nice axis labels.

%   Display the results :
%   max. horizontal distance & max. vertical height
% 

sprintf('%s','Maximum horizontal distance is ',max(x),' (m).')
sprintf('%s','Maximum vertical height is ',max(y),' (m).')

%
%   Plotting the result
%

figure(1); clf;                         % set up figure one, clear this figure if anything is on it
plot(x,y,'o');                          % plot trajectory -- y versus  x
xlabel('Horizontal Displ. (m)');        % axis labels must appear after the plot command
ylabel('Vertical Displ. (m)');          %axis limits must also appear after plot command
axis([0 x(iz1)*1.05 0 max(y)*1.05]);    % now limit the plot to those values of interest

The axis command limits the values of the plot to xmax to xmin. The value of x(iz1) is the maximum x trajectory.

Now this may have seemed long but really the entire code that you needed is only this:

%    inputs
%
g = 9.81;                    
N = 1000;                    
t_max = 10.0;                
v0 =     30.0                
theta = YOU CHOOSE!!  
%        
% calculations
%
t = 0:t_max/(N-1):t_max;      
x = YOU NEED TO PROVIDE THIS FORMULA; %  horizontal displacement
y = YOU NEED TO PROVIDE THIS FORMULA; %  vertical displacement
iz = find(y<0);
iz1 = iz(1);
%
%   Plotting the result
%
figure(1); clf;                         
plot(x(1:iz1),y(1:iz1),'o');            
xlabel('Horizontal Displ. (m)');        
ylabel('Vertical Displ. (m)');          
axis([0 x(iz1)*1.05 0 max(y)*1.05]);    

Only 13 lines of code.

7. The last thing you need to do is to learn how to loop through various values of theta to determine (in a brute force way) what the "best" angle is. The following code (with comments explaining what is going on) will do this for you. The two new Matalb concepts are looping ( FOR END loops and IF END logic) and multiple curve plotting. Create an "m-file"of the program listed below (for creating and running an "m-file" Click here). You will need to make the addition of the formulas for x and y dependence, otherwise this will loop through the theta values (you should vary the theta maximum and theta minumum values to satisfy yourself that you have the right region of initial angles).

%
%    Computer Assignment #1
%    Example problem 13-9, p49
%
%    Determine which angle of attack (30-60 deg.) gives the max. range
% 

clear all;  clf;

%
%    constants
%

g = 9.81;                     %  gravitational acceleration (m/s^2)
M = 7;                        %  No. of trials of angle of attack

%
%   Get the user inputs : Initial velocity & angle of attack
%

v0 = input('Initial velocity of the ball? (m/s) < 100 m/s  ');

theta0 = 30/180*pi;
thetaf = 60/180*pi;

j = 0;
for theta = theta0:(thetaf-theta0)/(M-1):thetaf  % loop through 7 values of theta
  j = j+1;
  
  t_max = 20;           %  final time of the flight; arbitray at this point.
  dt = 0.005;           % here we explicitly choose dt, the time increment.

  t = 0:dt:t_max;       %  time of flight 
  x = FORMULA;          %  horizontal displacement
  y = FORMULA;          %  vertical displacement

  iz = find(y<0);
  iz1 = iz(1);
  
  max_x(j) = x(iz1);  % array to find max value of x trajectory
  max_y(j) = max(y);  % max values of y trajectory

  figure(1); 
  
  if j == 1
    plot(x,y,'-');   hold on;   % hold on makes the plot stay up
  elseif j == 2
    plot(x,y,'--');  hold on;   % hold on allows for multiple curves  to be plotted
  elseif j == 3
    plot(x,y,'-.');  hold on;
  elseif j == 4
    plot(x,y,'r-');  hold on;
  elseif j == 5
    plot(x,y,'r-');  hold on;
  elseif j == 6
    plot(x,y,'r--'); hold on;
  elseif j == 7
    plot(x,y,'r-.'); hold on;
  end

end

axis([0 max(max_x)*1.05 0 max(max_y)*1.05]);  % plot to the maximum x value
h = legend('30 deg.','35 deg.','40 deg.','45 deg.','50 deg.',...
           '55 deg.','60 deg.');
ylabel('Vertical Displ. (m)');
xlabel('Horizontal Displ. (m)');