%
% hw1  sk  1/16/08
%
g = 9.8; %(m/s^2)
bom = 1;
II = 1;
K = 1;
xe(1) = pi/6; %(rad)
xe(2) = 0;
xe(3) = 10; % (rad/sec)
ue = K*cos(xe(1));
N = ceil(sqrt((g/cos(xe(1)))/xe(3)^2))
xe(3) = sqrt(g/cos(xe(1)))/N;
x0 = xe;

A = [ 0 1 0; 
     -g*(sin(xe(1)))^2/cos(xe(1)) -bom 2*g*sin(xe(1))/xe(3); 
     -K*sin(xe(1))/II 0 0]
  disp('Eigenvalues of A')
eig(A)
B = [0 0 -1/II]'

  t = [0:.01:10]';

npts = max(size(t));
sys = ss(A,B,eye(3),0);
y = lsim(sys,ue*ones(size(t)),t,x0);

figure(1)
  plot(t,y(:,1)+x0(1)*ones(size(t)))
  ylabel('\phi_0 + \Delta\phi (rad)')
  xlabel('Time (s)')
  title('Centrifugal Governor Vertical Angle ') 
  grid
figure(11)
  plot(t,y(:,1))
  ylabel('\Delta\phi (rad)')
  xlabel('Time (s)')
  grid

plot(t,y(:,2))
figure(3)
plot(t,y(:,3))
grid
