%
%   HW9_DEMO3.M  Compute motorboat characteristics versus time (HW9_DEMO3)
%
%   This problem was derived from Probs. 2&3 pg. 103 in the 2nd Ed. of Edwards and 
%   Penny (Prentice Hall, 2000).  It addresses the NUMERICAL solution of a 
%   particular dynamics problem that models the speed and position of a motorboat 
%   as it coasts to a stop.
%
%   The force balance equation for this problem is:    dv/dt = -k*v    
%   and the position equation is:                      dx/dt = v   
%   where k is a constant that needs to be determined based on data for the 
%   particular problem.  We will solve the problem for a range of k's.
%
%   Calls HW9_DEMO3_EQN.M to evaluate the desired vector function, f(z,t).  Here z is
%   a 2-element vector with z(1) = x and z(2) = v.  With this definition, the 
%   above equations can be written as
%       dz(1)/dt = z(2)
%       dz(2)/dt = -k*z(2)
%
%   File prepared by J. R. White, UMass-Lowell (March 1999)
%    --> modfied slightly for 2nd Ed.  (Oct. 1999)
%    --> changed names for consistency (Jan. 2000)
%    --> minor changes (March 2002)
%
      
%
%   getting started
      clear all;   close all;   nfig = 0;
      global kk              % constant to be used in ODE function file
%
%   define integration variables and initial condition
      to = 0;   tfinal = 150;% time in seconds
      t = linspace(to,tfinal,151);% this specifies specific time points for the ode soln
      xo = 0;   vo = 40;     % initial position (ft) and speed (ft/s)
      zo = [xo vo]';         % initial conditions in vector form
%
%   evaluate x(t) and v(t) (ie. components of vector z) for various values of kk
      cols = ['r- ';'g: ';'b-.';'m--'];
      n = 0;  kk = input('Input an initial guess for k (proportionality constant): ');
      while kk > 0
        n = n+1;   if n > 4,  n = 1;   end
        [t,z] = ode23('hw9_demo3_eqn',t,zo);
        subplot(2,1,1),plot(t,z(:,1),cols(n,:),'LineWidth',2),grid on,hold on
        nn = min(find(t > 80));  xt = 122;  yt = z(nn,1)-40;
        text(xt,yt,['k = ',num2str(kk)])
        subplot(2,1,2),plot(t,z(:,2),cols(n,:),'LineWidth',2),grid on,hold on
        subplot(2,1,2),plot(10,20,'bx','LineWidth',2),grid on,text(10.5,21,'expt') % mark expt pt
        kk = input('Input another value for k (enter zero to quit): ');
      end
%   label subplot #1
      subplot(2,1,1),ylabel('Position (ft.)')
      title('Motorboat Dynamics while Coasting to a Stop (HW9 DEMO3)')
      hold off
%   label subplot #2
      subplot(2,1,2),xlabel('Time (sec)'),ylabel('Speed (ft./sec)')
      hold off
%
%   end of problem
%
