% Script for ME 380 Problem 19(a) with the data I gave in the problem.
% This performs the state-space modeling of an armature voltage-controlled
% DC motor.  The state vector is x = [i_a; theta_m; omega_m], and the input
% is u = e_a.  The output is y = [i_a; theta_m (rev); omega_m (rpm)].  The
% input is a voltage pulse.

% 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)

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

A = [-R_a/L_a  0  -K_b/L_a;     % 3x3 A system matrix
         0     0      1;
      K_t/J_m  0  -b_m/J_m];
  
B = [1/L_a; 0; 0];  % 3x1 B input matrix

C = [1     0         0;         % Output matrix C produces i_a (A),
     0  1/(2*pi)     0          % theta_m (revolutions), and
     0     0     60/(2*pi)];    % omega_m (rpm)
 
D = [0;0;0];    % Feed-thru matrix D is zero for this physical system

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 = lsim(motor,u,t);    % Perform simulation

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)