%
%   PROJECTILE2D_1.M   Plot object trajectories for simple 2-D motion
%                     assuming only a gravity force in the y-direction
%
%   Projectile motion is governed by Newton's 2nd Law.  If the only force present
%   is the downward force of gravity (i.e. no air friction), then the resulting 
%   equations of motion are quite simple.  In particular, let's assume an x-y 
%   coordinate system as shown below:
%         ^
%      y  |
%         |            where +y is up and +x is to the right
%         ------>
%              x
%
%   We will assume that the object (a bottle cap, for example) has an initial 
%   position given by  
%          x(0) = xo = 0    and     y(0) = yo  (user specified)
%   
%   The initial velocity is given by the object's speed, Vo, and its angle, theta0,
%   relative to the horizontal axis.  With these variables specified, the initial
%   x-directed velocity, vxo,  and initial y-directed velocity, vyo, are given by
%          vxo = Vo*cos(theta0)    and    vyo = Vo*sin(theta0)
%
%   For the case of no air friction, Newton's law gives
%      x position vs time:    x(t) = vxo*t
%      y position vs time:    y(t) = yo + vyo*t - (g/2)*t^2
%
%   Thus, given yo, Vo, and theta0, we can easily evaluate these expressions for several
%   values of t and plot x(t) vs t and y(t) vs t -- or, of even more interest, is a plot 
%   of y(t) vs x(t), since this shows the projectile's trajectory in the x-y plane.  
%
%   Now to make this interesting, let's say you are setting at your desk and you open 
%   a bottle of soda or beer (or whatever).  However, your garbage can is across the room.
%   You decide to "give it a try" anyway, and attempt to toss the bottle cap into the 
%   trash can.  The conditions are as follows:
%    -- vertical position of cap with you sitting with arm extended is about 4 ft or 1.3 m
%    -- location of front edge of can is about 10 ft or 3 m from where you are sitting
%    -- the can is about 1.5 ft or 0.45 m high with a depth of about 1 ft or 0.3 m
%    -- a reasonable initial speed for tossing the cap is about 20 - 30 ft/s or 6 - 10 m/s
%    -- the initial angle has to be such that the x-directed velocity is positive which
%       gives a range of -pi/2 < theta0 < pi/2. This is the maximum range possible and our
%       goal will be to determine the values of theta0 that lead to success.
%   All these values can vary somewhat, but they represent a reasonable set of conditions.
%
%   Now, with a set of specific conditions chosen as indicated above, the goal here is to 
%   find the range of angles which lead to you successfully getting the bottle cap into  
%   the garbage can.  Good Luck...
%
%   File prepared by J. R. White, UMass-Lowell (last update: September 2017)
%

      clear all;   close all;   nfig = 0;
%
%   define some initial parameters (these can vary somewhat for different cases)
      g = 9.8;        % gravitational acceleration constant (m/s^2)
      x1 = 3;         % x location of front edge of can (m)
      x2 = x1+0.3;    % x location of back edge of can (m)
      y1 = 0.45;      % y location of top edge of can (m)
      yo = 1.3;       % initial height of bottle cap (m)
      Vo = 6;         % initial speed of bottle cap (m/s)
%
%   define range of initial angles (in radians)  -- note that we can modify the low 
%   end of the theta0 range from simple geometry considerations
      thmin = atan((yo-y1)/x1);   % use this as the min theta instead of -pi/2
      Nth = 50;  theta0 = linspace(-thmin,pi/2,Nth);
%
%   define time domain for simulation (it should only take a few seconds for the cap 
%   to hit the floor --  we will simply ignore any negative y values)
      Nt = 41;   t = linspace(0,2,Nt)';
%
%   now evaluate desired quantities 
%
%   Note: "dot arithmetic" is used here to evaluate the desired quantities at Nt 
%   different times in a single statement -- this vector processing capability
%   is one of the many very nice features available in Matlab.  However, we also  
%   want to do all these different calculations for Nth different values of theta0.
%   This will be handled with a looping structure.  When finished, we will have a 
%   matrix of x and y values --> since they are functions of two variables, t = time
%   and theta0 = initial angle of bottle cap.
%
      x = zeros(Nt,Nth);   y = zeros(Nt,Nth);     % initialize variables
      for i = 1:Nth                               % do the calculations
        th = theta0(i);
        vxo = Vo*cos(th);   vyo = Vo*sin(th);
        x(:,i) = vxo*t;
        y(:,i) = yo + vyo*t -(g/2)*t.*t;
      end    
%
%   plot y(t) vs x(t) (for all values of thetha0) along with an outline of the can
      nfig = nfig+1;   figure(nfig)
      plot(x,y), grid, hold on
      plot([x1 x1],[0 y1],'k-',[x2 x2],[0 y1],'k-',[x1 x2],[y1 y1],'k-','LineWidth',2)
      axis([0 4 0 4])
      title(['Projectile2d\_1:  Projectile Trajectories for many \theta_o Values (V_o = ', ...
             num2str(Vo),' m/s)'])
      xlabel('x location (m)'),ylabel('y location (m)')
      hold off
%
%  from above, we can see that there are two ranges of angles that lead to success
%    direct path -- low angle trajectory that leads to the shortest path to the can
%    indirect path -- high angle path that takes a longer time to reach the can
%
%  for the current case with Vo = 6 m/s, it appears that the low angle paths are at  
%  about 5 - 15 degrees and the high angle successful trajectories are in the 60 - 70  
%  degree range.
%
%  to explore these specific angles, let's put in an interactive loop and ask the user 
%  for specific values of theta0
%
%  first let's draw the x-y grid and the garbage can in a new figure window
      nfig = nfig+1;   figure(nfig)
      plot(0,yo,'ks'), grid, hold on
      plot([x1 x1],[0 y1],'k-',[x2 x2],[0 y1],'k-',[x1 x2],[y1 y1],'k-','LineWidth',2)
      axis([0 4 0 4])
      title(['Projectile2d\_1:  Projectile Trajectories in X-Y Plane (V_o = ', ...
             num2str(Vo),' m/s)'])
      xlabel('x location (m)'),ylabel('y location (m)')
%
%   now loop over user-specified values of theta0 and plot x-y trajectory for each
      cont = 1;
      while cont == 1
        th_d = input('Input a value for theta0 in degrees: ');
        th = th_d*pi/180;
        vxo = Vo*cos(th);   vyo = Vo*sin(th);
        xx = vxo*t;  yy = yo + vyo*t -(g/2)*t.*t;
        plot(xx,yy,'LineWidth',2),gtext(['\theta_o = ',num2str(th_d),' ^o'])
        cont = menu('Perform another simulation?', ...
                    'Yes -- I want to try another angle', ...
                    'No Thanks -- I have seen enough');
      end
      hold off
%
%   end of program
