% Script for ME 380 Problem 19(c) with the data I gave in the problem.
% This is a simplified motor model with armature inductance neglected.
% The state vector is now  x = [theta_L; omega_L], but we are still
% intereted in output y = [i_a; theta_L; omega_L].  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)

rho = (2.7/1000)*(100)^3;   % Density of Aluminum (kg/m^3)
r_L = 100/1000;	% Radius of load disk (m)
t_L = 10/1000; 	% Thickness of load disk (m)
b_L = 4/1000;   % Load damping (N-m-s/rad)

n = 50;     % Gear ratio (speed reducer)

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

% Find mass moment of inertia of load disk.

V_L = (pi*r_L^2)*t_L;	% Volume of load disk (m^3);
m_L = rho*V_L;          % Mass of load disk (kg);
J_L = (1/2)*m_L*r_L^2;  % Mass moment of inertia of load disk

% Find the total inertia and damping as felt by the motor

J_t = J_m+J_L/n^2;  % Total inertia (kg-m^2)
b_t = b_eq+b_L/n^2; % Total damping (N-m-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_t/J_t];
  
B = [0; K_t/(J_t*R_a)];  % 2x1 B input matrix


C = [   0      -K_b/R_a;        % Output matrix computes current i_a, and
    180/(n*pi)     0;           % changes units for theta_L (deg), and
        0      180/(n*pi)];     % omega_L (deg/s)
 
D = [1/R_a;0;0];    % Feed-thru matrix D computes part of 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 motor current (computed)
figure;
plot(t,y(:,2));     % Plot first output (theta_L in deg)
figure;
plot(t,y(:,3));     % Plot second output (omega_L in deg/s)