% Script for ME 380 Problem 19(b) with the data I gave in the problem.
% This is a simplified motor model with armature inductance neglected.
% The state vector is now just x = [theta_m; omega_m], but we are still
% intereted in output y = [i_a; theta_m; omega_m].  We can be clever and
% form current i_a during the simulation by proper definition of the C and
% D matrices.

% Motor parameter values

J_m = 1.54e-6;  % Motor armature inertia (kg-m^2)
b_m = 2.7e-6;   % Motor damping (N-m-s/rad)
L_a = 0.69e-3;  % Motor armature inductance (H)
R_a = 1.53;     % Motor armature resistance (ohms)
K_t = 22.3e-3;  % Motor torque constant (N-m/A)
K_b = 0.0223;   % Motor back EMF constant (V-s/rad)

b_eq = ((K_t*K_b)/R_a)+b_m;   % "Equivalent" damping coefficient

% Next define the state-space "A B C D" matrices

A = [0      1;     % 2x2 A system matrix
     0  -b_eq/J_m];
  
B = [0; K_t/(J_m*R_a)];  % 2x1 B input matrix

% Now...in addition to changing units for theta_m and omega_m, the C and
% D matrices can be used to calculate motor current i_a.  Since we really
% have 3 outputs, 2 state variables, and 1 input, matrix C is of dimension
% 3x2, and matrix D is of dimension 3x1.

C = [   0     -K_b/R_a;     % Output matrix computes i_a (A), and changes
     1/(2*pi)     0         % units of theta_m to (rev) and
        0     60/(2*pi)];    % omega_m to (rpm)
 
D = [1/R_a; 0; 0];    % Feed-thru matrix D used to compute current i_a

motor = ss(A,B,C,D);    % Construct the state-space model of the motor

% Next construct the input u, which consists of a 10V pulse of duration 0.1
% second, with a total simulation time of 0.2 seconds.  As suggested, use a
% time step of dt = 0.001 second (1 msec).
%
% The input will consist of two parts: u1, which will be the 10V portion;
% it will be 101 samples in length, and u2, which will be the "zero"
% portion; it will be 100 samples in length.  These will be formed as
% column vectors, and "stacked" together.

u1 = ones(101,1)*10;    % 101 samples of magnitude 10
u2 = zeros(100,1);      % 100 samples of magnitude zero

u = [u1;u2];    % Stack them on "top" of each other

% The next step is to simulate the motor to the pulse input.  The syntax of
% the "lsim" function is: y = lsim(sys,u,t).  So we need a time vector t
% corresponsponding to the input u.

dt = 0.001;     % Time step 1 msec
t = [0:dt:0.2]';  % Time vector from 0 to 0.2 sec (transpose to get column)

% Finally, perform the simulation and plot the results.

[y,t,x] = lsim(motor,u,t);    % Perform simulation; get state vector

plot(t,y(:,1));     % Plot first output (i_a in A)
figure;
plot(t,y(:,2));     % Plot second output (theta_m in rev)
figure;
plot(t,y(:,3));     % Plot third output (omega_m in rpm)