% VehSusp_script.m -->  with VehSuspBlock.mdl
% EGR 390 HW 4
% J Cardell, February 2009
% Script for Simulink model of the vehicle suspension system

% Constants
m1 = 100;   % mass of tire
k1 = 50;    % spring constant for tire assembly
m2 = 4000;  % mass of vehicle
k2 = 75;    % spring constant for suspension system
b = 50;     % damping coefficient shock absorber

% A, B, C, D system matrices
% x1 = omega, angular velocity;  x2 = Tk, spring torque
% The system starts with stored energy in the mass and spring
A = [-b/m1 1/m1 b/m1 -1/m1; -k1 0 0 0 ; b/m2 0 -b/m2 1/m2; k2 0 -k2 0];  
B = [1/m1 ; 0; 0; 0];
C = [0 -1/k1 0 0; 0 -1/k1 0 -1/k2];  % Vertical position selected as output
D = [0];

% Initial conditions could be included, via the state space block in 
% Simulink, but since there is an input force, via the uneven road surface,
% we can observe the system behavior without included initial stored energy

% Run the simulink model
sim('HW4_2_VehSuspBlock');

% Define data vectors written to Matlab workspace by Simulink 'To
% Workspace' block.  Note that the 'C' matrix identified 2 output
% variables, the veritical displacement of the tire and of the bus
t = road.time;  % The 'time' vector, used for the x-axis values
u = road.signals.values; % Input data vector
y1 = vert_position.signals.values(:,1);  % Output - vertical position
y2 = vert_position.signals.values(:,2);

% plot the data
figure
hold on

subplot(3, 1, 1)
plot(t, u, 'c-', 'LineWidth', 2)
title('Input road surface')

subplot(3, 1, 2)
plot(t, y1, 'r-', 'LineWidth', 2)
title('Response of a Vehicle Suspension System to an Uneven Road Surface')
title('Tire response: Vertical displacement')
ylabel('Displacement (m)')

subplot(3, 1, 3)
plot(t, y2, 'b-', 'LineWidth', 2)
title('Vehicle response: Vertical displacement')
xlabel('t (sec)')
ylabel('Displacement (m)')
grid
hold off