File: ef230-2023-08/www/examples/simulation_hibbler_12_19.m Download
%! Simulation of a car 'crash' with extensive plotting
% University of Tennessee : Chris Pionke
%========================================================================================%
%                          |                          |                          |       %
% EF 102                   |           2/28/04        |  Merlot Red Z4           | Pg 1  %
%                          |                          |                          |       %
%========================================================================================%
%                                                                                        %
% This program is a simulation of problem 12-19 in "Dynamics" 9th ed by R.C Hibbeler.    %
% Car A is replaced by BMW Z4 and Car B is replaced by a Mini Cooper.                    %
%                                                                                        %
% The user has the ability to input (i.e., change) the values of various parameters in   %
% problem. The parameters that the user can specify are:                                 %
%         the initial separation distance                                                %
%         the reaction time of driver B                                                  %
%         the intimal velocity of the two cars (they are the same)                       %
%         the Z4 deceleration                                                            %
%         the Mini deceleration                                                          %
%                                                                                        %
% All output is sent to figure 1.                                                        %
% All variable names are defined in the in-line comments.                                %
% The files "z4.jpg" and "mini.jpg" must be in the current directory.              %
%                                                                                        %
%========================================================================================%
% Clear memory and the Command Window.
clear all; clc;

% Set default values for the user input. They are the initial separation distance (gap),
% the reaction time (rt), the intimal velocity of the cars (v_ini), the Z4 deceleration
% (a_Z4), and the Mini deceleration (a_mini).
gap    = 15;
rt     = 0.75;
v_ini  = 60;
a_Z4   = 12;
a_mini = 15;

% Get user input with a dialog box and check for errors.
check = 0;
while check == 0
    values = {num2str(gap),num2str(rt),num2str(v_ini),num2str(a_Z4),num2str(a_mini)};
    values = inputdlg(...
        {'Initial distance between cars - ft (must be >= 1 ft)', ...
         'Reaction time - s (must be >= 0.1 s)',...
         'Initial speed of both cars - ft/s  (must be <= 100 ft/s)', ...
         'Z4 deceleration - ft/s^2  (must be <= 20 ft.s^2)', ...
         'Mini deceleration - ft/s^2  (must be <= 20 ft.s^2)'}, ...
         'Hibbeler Problem 12-19',1,values);
    if length(values)==5
        check  = 1;
        gap    = str2num(char(values(1)));
        rt     = str2num(char(values(2)));
        v_ini  = str2num(char(values(3)));
        a_Z4   = str2num(char(values(4)));
        a_min  = str2num(char(values(5)));

        if gap < 1 
            waitfor(errordlg('Invalid Gap Distance'));
            check = 0;
        end
        if rt < 0.1 
            waitfor(errordlg('Invalid Reaction Time'));
            check = 0;
        end
        if v_ini < 0 | v_ini > 100
            waitfor(errordlg('Invalid Initial Speeds'));
            check = 0;
        end
        if a_Z4 < 0 | a_Z4 > 20
            waitfor(errordlg('Invalid Z4 Deceleration'));
            check = 0;
        end
         if a_min < 0 | a_mini > 20
            waitfor(errordlg('Invalid Mini Deceleration'));
            check = 0;
        end
       if a_mini <= a_Z4
            waitfor(errordlg('Mini Decleeration Must Exceed Z4 Speed'));
            check = 0;
        end
    end
end

%========================================================================================%

% Set initial speeds, initial positions, decellerations, weights and masses of both cars.
v_A  = v_ini;
v_B  = v_ini;
a_A  = a_Z4;
a_B  = a_mini;
x_B  = 0;
x_A  = x_B + gap;
wt_A = 3000;
wt_B = 2500;
g    = 32.2;
m_A  = wt_A ./ g;
m_B  = wt_B ./ g;

% Calculate the stopping time and position of both cars and determine the maximum stopping
% time and position.
tstop_A = v_A ./ a_A;
tstop_B = (v_B ./ a_B) + 0.75;
tmax = max([tstop_A tstop_B]);
sstop_A = x_A + v_A .* tstop_A - (a_A ./ 2) .* tstop_A.^2;
sstop_B = x_B + v_B .* tstop_B - (a_B ./ 2) .* tstop_B.^2;
smax = max([sstop_A sstop_B]);

% Set the simulation time increment and create all the time values.
dt = 0.01;
t = 0:dt:tmax;

% Set the initial values for the position vectors (xA & xB), the velocity vectors
% (vA & vB), and the kinetic energy vectors (keA, keB, & keT). 
xA(1) = x_A;
xB(1) = x_B;
vA(1) = v_A;
vB(1) = v_B;
keA(1) = 0.5 .* m_A .* vA(1).^2;
keB(1) = 0.5 .* m_B .* vB(1).^2;
keT(1) = keA(1) + keB(1);

%========================================================================================%

% Activate figure 1, clear all parameters and set the figure's position and size.
figure(1);
clf reset;
set(1,'Position',[100 100 800 560]);
% Note, if you are using a 800x600 resolution screen, use the following (remove the %)
% set(1,'Position',[0 30 800 560]);

% Plot the background rectangles for all the "meters" and plots.
axes('Position',[.04 .82 .225 .2],'Color',[1 .82 .4]);
hold on;
box on;
set(gca,'YTick',[ ],'XTick',[ ])

axes('Position',[.04 .705 .225 .105],'Color',[1 .82 .4]);
hold on;
box on;
set(gca,'YTick',[ ],'XTick',[ ])

axes('Position',[.04 .56 .225 .135],'Color',[1 .82 .4]);
hold on;
box on;
set(gca,'YTick',[ ],'XTick',[ ])

axes('Position',[.285 .56 .33 .44],'Color',[1 .82 .4]);
hold on;
box on;
set(gca,'YTick',[ ],'XTick',[ ])

axes('Position',[.635 .56 .33 .44],'Color',[1 .82 .4]);
hold on;
box on;
set(gca,'YTick',[ ],'XTick',[ ])

axes('Position',[.04 .05 .38 .365],'Color',[.7 .7 .5]);
hold on;
box on;
set(gca,'YTick',[ ],'XTick',[ ])

axes('Position',[.43 .05 .535 .15],'Color',[.7 .7 .5]);
hold on;
box on;
set(gca,'YTick',[ ],'XTick',[ ])

axes('Position',[.43 .22 .26 .195],'Color',[.7 .7 .5]);
hold on;
box on;
set(gca,'YTick',[ ],'XTick',[ ])

axes('Position',[.7 .22 .265 .195],'Color',[.7 .7 .5]);
hold on;
box on;
set(gca,'YTick',[ ],'XTick',[ ])

%========================================================================================%

% Set up the various "meters" and "plots" that are above the vehical simulation window.

% Set up the "Position Meter."
axes('Position',[.05 .835 .2 .12]);
hold on;
box on;
axis([0 1 0 1])
set(gca,'YTick',[ ],'XTick',[ ])
title('Position Meter');
text(.05,.8,['Mini  -        ft'],'Color','r','FontName','FixedWidth',...
                   'EraseMode','xor');
text(.05,.5,['Z4    -        ft'],'Color','b','FontName','FixedWidth',...
                   'EraseMode','xor');
text(.05,.2,['Apart -        ft'],'Color','k','FontName','FixedWidth',...
                   'EraseMode','xor');
tpB = text(.78,.8,[num2str(xB(1),'%6.2f')],'Color','r','FontName','FixedWidth',...
                   'HorizontalAlignment','right','EraseMode','xor');
tpA = text(.78,.5,[num2str(xA(1),'%6.2f')],'Color','b','FontName','FixedWidth',...
                   'HorizontalAlignment','right','EraseMode','xor');
tpT = text(.78,.2,[num2str(xA(1)-xB(1),'%6.2f')],'Color','k','FontName','FixedWidth',...
                   'HorizontalAlignment','right','EraseMode','xor');

% Set up the "Time Meter."
axes('Position',[.05 .715 .2 .05]);
hold on;
box on;
axis([0 1 0 1]);
set(gca,'YTick',[ ],'XTick',[ ]);
title('Time Meter');
text(.15,.4,['t =       sec'],'FontName','FixedWidth','EraseMode','xor');
tt = text(.56,.4,[num2str(t(1),'%5.2f')],'FontName','FixedWidth',...
           'HorizontalAlignment','right','EraseMode','xor');
              
% Set up the "Speed Meter."
axes('Position',[.05 .57 .2 .08]);
hold on;
box on;
axis([0 1 0 1])
set(gca,'YTick',[ ],'XTick',[ ])
title('Speed Meter');
text(.05,.7,['Mini -      ft/sec'],'Color','r','FontName','FixedWidth','EraseMode','xor');
text(.05,.3,['Z4   -      ft/sec'],'Color','b','FontName','FixedWidth','EraseMode','xor');
tvB = text(.63,.7,[num2str(vB(1),'%4.1f')],'Color','r','FontName','FixedWidth',...
                   'HorizontalAlignment','right','EraseMode','xor');
tvA = text(.63,.3,[num2str(vA(1),'%4.1f')],'Color','b','FontName','FixedWidth',...
                   'HorizontalAlignment','right','EraseMode','xor');

% Set up the "Position Plot."
axes('Position',[.35 .65 .25 .3]);
hold on;
box on;
ppB = plot(t(1),xB(1),'r','EraseMode','none');
ppA = plot(t(1),xA(1),'b','EraseMode','none');
axis([0 tmax 0 smax.*1.1]);
title('Position Plot');
xlabel('Time - sec');
ylabel('Position - ft');
legend('Mini','Z4',4)

% Set up the "Speed Plot."
axes('Position',[.7 .65 .25 .3]);
hold on;
box on;
pvB = plot(t(1),v_B,'r','EraseMode','none');
pvA = plot(t(1),v_A,'b','EraseMode','none');
axis([0 tmax 0 v_ini+10]);
title('Speed Plot');
xlabel('Time - sec');
ylabel('Speed - ft/sec');
legend('Mini','Z4',1);

%========================================================================================%

% Set up the various "meters" and "plots" that are below the vehical simulation window.

% Set up the "Kinetic Energy Plot."
axes('Position',[.1 .12 .3 .25]);
box on;
hold on;
pkeB = plot(t(1),keB(1),'r','EraseMode','none');
pkeA = plot(t(1),keA(1),'b','EraseMode','none');
pkeT = plot(t(1),keT(1),'k','EraseMode','none');
axis([0 tmax 0 keT(1).*1.1]);
title('       Kinetic Energy Plot');
xlabel('Time - sec');
ylabel('Energy - ft\cdotlb');
legend('Mini','Z4','Total',-1);

% Set up the "Kinetic energy Meter."
axes('Position',[.45 .25 .22 .12]);
hold on;
box on;
axis([0 1 0 1])
set(gca,'YTick',[ ],'XTick',[ ])
title('Kinetic Energy Meter');
text(.05,.8,['Mini  -         ft\cdotlb'],'Color','r','FontName','FixedWidth',...
                'EraseMode','xor');
text(.05,.5,['Z4    -         ft\cdotlb'],'Color','b','FontName','FixedWidth',...
                'EraseMode','xor');
text(.05,.2,['Total -         ft\cdotlb'],'Color','k','FontName','FixedWidth',...
                'EraseMode','xor');
tkeB = text(.74,.8,[num2str(keB(1),'%6.0f')],'Color','r','FontName','FixedWidth',...
                'HorizontalAlignment','right','EraseMode','xor');
tkeA = text(.74,.5,[num2str(keA(1),'%6.0f')],'Color','b','FontName','FixedWidth',...
                'HorizontalAlignment','right','EraseMode','xor');
tkeT = text(.74,.2,[num2str(keT(1),'%6.0f')],'Color','k','FontName','FixedWidth',...
                'HorizontalAlignment','right','EraseMode','xor');

% Set up the "Initial Conditions Box."
axes('Position',[.71 .23 .25 .14]);
hold on;
box on;
axis([0 1 0 1]);
set(gca,'YTick',[ ],'XTick',[ ]);
title('Initial Conditions');
text(.05,.8,['Gap distance = ',num2str(gap,'%4.0f'),' ft'],'FontName','FixedWidth');
text(.05,.5,['a_{Mini} = -',num2str(a_mini,'%4.0f'),' ft/s^2'],'FontName','FixedWidth',...
           'Color','r');
text(.05,.2,['a_{Z4  } = -',num2str(a_Z4,'%4.0f'),' ft/s^2'],'FontName','FixedWidth',...
           'Color','b');
       
% Set up the "Assumptions Box."
axes('Position',[.45 .07 .5 .07]);
hold on;
box on;
axis([0 1 0 1]);
set(gca,'YTick',[ ],'XTick',[ ]);
title('Assumptions');
text(.02,.25,['The simulation will stop if the cars collide.'],'FontName','FixedWidth');
text(.02,.65,['Weight_{Mini} = ',num2str(wt_B,'%04.0f'),' lb'],'FontName','FixedWidth',...
           'Color','r');
text(.56,.65,['Weight_{Z4} = ',num2str(wt_A,'%04.0f'),' lb'],'FontName','FixedWidth',...
           'Color','b');

%========================================================================================%

% Set up the "Vehical Simulation Window."
axes('Position',[.04 .45 .925 .1]);
hold on;
box on;
axis([-50 smax*1.25 0 10])

% Add the text "ft" to the tick labels.
hhh = get(gca,'XTickLabel');
ggg = size(hhh);
for iii = 1:1:ggg(1)
    hhh(iii,ggg(2)+1) = ' ';
    hhh(iii,ggg(2)+2) = 'f';
    hhh(iii,ggg(2)+3) = 't';
end
set(gca,'YTick',[ ])
set(gca,'XTickLabel',hhh);

% Read in the figure files. These two files must be in the current directory.
urlwrite('http://ef.engr.utk.edu/ef230-2009-08/modules/misc/z4.jpg','z4.jpg');
urlwrite('http://ef.engr.utk.edu/ef230-2009-08/modules/misc/mini.jpg','mini.jpg');
z4img = imread('z4.jpg','jpg');
miniimg = imread('mini.jpg','jpg');

% Plot the pictures of the two cars in their initial positions.
image([xA(1) xA(1)+0.13333.*(1.25.*smax+50)],[0 8],flipdim(z4img,1),'EraseMode','xor');
image([xB(1)-0.11111.*(1.25.*smax+50) xB(1)],[0 8],flipdim(miniimg,1),'EraseMode','xor');


%========================================================================================%
button = questdlg('Run Simulation?',...
'Continue','Run','Run');
% Loop over all time values of the simulation starting from the 2nd value.
for ii = 2:1:length(t)
    
    % For the current time step, calculate the new postion of both cars based on their
    % position at the previous time step and their velocity during the current time step.
    xA(ii) = xA(ii-1) + vA(ii-1) .* dt;
    xB(ii) = xB(ii-1) + vB(ii-1) .* dt;
    
    % Calculate the velocity of car A at the end of the current time step. The new value
    % is based on the velocity at the beginng of the time step and its deceleration.
    vA(ii) = vA(ii-1) - a_A .* dt;
    
    % Calculate the velocity of car B at the end of the current time step. The new value
    % is based on the velocity at the beginng of the time step and its deceleration.
    % If time is <= the "reaction time," the velocity of B is constant and is set to
    % its initial velocity.
    if t(ii) <= rt
        vB(ii) = vB(ii-1);
    else
        vB(ii) = vB(ii-1) - a_B .* dt;
    end
    
    % Check both cars to see if they have come to a stop. If so, set their velocity's
    % to zero.
    if vA(ii) < 0
        vA(ii) = 0;
    end
    if vB(ii) < 0
        vB(ii) = 0;
    end
    
    % Calculate the kinetic energies at the current time step.
    keA(ii) = 0.5 .* m_A .* vA(ii).^2;
    keB(ii) = 0.5 .* m_B .* vB(ii).^2;
    keT(ii) = keA(ii) + keB(ii);
    
    % Check to see if the car B "crashes" into car A. If so, ouput a message to the
    % figure and stop the simulation.
    if xB(ii) >= xA(ii)
        text(xB(ii),11,['CRASH'],'BackgroundColor','w','Color','r','FontSize',24,...
            'FontWeight','bold','HorizontalAlignment','center','EraseMode','xor');
        drawnow
        break
    end
    
    % If the simulation has not stopped, update all "Plots, Meters," and the
    % "Vehicle Simulation Window."
    set(ppB,'XData',[t(ii-1) t(ii)],'YData',[xB(ii-1) xB(ii)])
    set(ppA,'XData',[t(ii-1) t(ii)],'YData',[xA(ii-1) xA(ii)])
    set(pvB,'XData',[t(ii-1) t(ii)],'YData',[vB(ii-1) vB(ii)])
    set(pvA,'XData',[t(ii-1) t(ii)],'YData',[vA(ii-1) vA(ii)])
    set(pkeB,'XData',[t(ii-1) t(ii)],'YData',[keB(ii-1) keB(ii)])
    set(pkeA,'XData',[t(ii-1) t(ii)],'YData',[keA(ii-1) keA(ii)])
    set(pkeT,'XData',[t(ii-1) t(ii)],'YData',[keT(ii-1) keT(ii)])
    set(tkeB,'String',[num2str(keB(ii),'%6.0f')]);
    set(tkeA,'String',[num2str(keA(ii),'%6.0f')]);
    set(tkeT,'String',[num2str(keT(ii),'%6.0f')]);
    set(tpB,'String',[num2str(xB(ii),'%5.2f')]);
    set(tpA,'String',[num2str(xA(ii),'%5.2f')]);
    set(tpT,'String',[num2str(xA(ii)-xB(ii),'%5.2f')]);
    set(tvB,'String',[num2str(vB(ii),'%4.1f')]);
    set(tvA,'String',[num2str(vA(ii),'%4.1f')]);
    set(tt,'String',[num2str(t(ii),'%5.2f')]);
    image([xA(ii) xA(ii)+0.13333.*(1.25.*smax+50)],[0 8],flipdim(z4img,1),...
           'EraseMode','xor');
    image([xB(ii)-0.11111.*(1.25.*smax+50) xB(1) xB(ii)],[0 8],flipdim(miniimg,1),...
           'EraseMode','xor');
    drawnow
end

% End of program.