Problem Statement: Using Matlab, plot the x-y trajectory of a mass flying through the air. Assume that only a gravity force acts on the mass, and that the initial velocity and initial angle (measured from the horizontal) are given. The initial position is x=0, y=h. By plotting trajectories, determine the following:

1. Use an initial velocity of v0=30 m/s and an initial height of h=20 m. On a single figure, plot the x-y trajectory for various initial angles: 30, 35, 40, 45, 50, 55, and 60 degrees. From this figure, determine the (approximate) angle which yields the longest distance in the horizontal plane (x-direction distance).
2. Does this best angle depend on the initial height? On another figure, plot the x-y trajectories for the same initial velocity and the same initial angles as in part 1. However, use an initial height of h=0 (i.e., the projectile starts from the ground). What is the best angle for this case? Does this result agree with the solution for homework problem 2.68?
3. Now use the Matlab plot commands to plot the x and y dependence on time. Plot x vs. t and y vs. t for a 30 degree initial angle with v0 = 30 m/s and h = 20 m.
4. Print out your "m-files" that perform the tasks above.

• First you should familiarize yourself with the flight of a projectile (mass) which has been given an initial velocity and direction with only gravity acting on the mass. All formulas you need were given in class.
• Follow the Matlab Tutorial below for instructions on how to complete the above assignment.  

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 v0 = 30 m/s.
v0 = 30.0;              % Initial velocity of the ball (m/s)
h = 20.0;               % Initial height of the ball (m/s) --> you can change
theta = YOU CHOOSE!! ;  % initial angle in RADIANS = pi*(degrees)/180

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 provide 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, y 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   % DO NOT 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.
4. Now that the x and y dependence is known, we can plot the variables.

plot(x,y,'-');

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 an 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),'-');

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).
6. Now you can plot the results. The following prints a few pertinent 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,'-');                          % 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;                
h = 20.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),'-');            
xlabel('Horizontal Displ. (m)');        
ylabel('Vertical Displ. (m)');          
axis([0 x(iz1)*1.05 0 max(y)*1.05]);    

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 Matlab 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, but this will loop through the theta values for you.

%
%    Computer Assignment #1
%    ME 240
%
%    Determine which initial angle (30-60 deg.) gives the max. range
% 

clear all;  clf;

%
%    constants
%

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

%
%   Get the user input : Initial height
%

h = input('Initial height of the ball? (m) ');

%
% loop through 7 values of theta
%

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

j = 0;
for theta = theta0:(thetaf-theta0)/(M-1):thetaf  
  j = j+1;
  
  t_max = 10;           %  final time of the flight; arbitrary 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
legend('30 deg.','35 deg.','40 deg.','45 deg.','50 deg.',...
           '55 deg.','60 deg.');
ylabel('Vertical Displ. (m)');
xlabel('Horizontal Displ. (m)');